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A  time-dependent  molecular  orbital  method  has  been  developed  to  study 
charge  transfer  in  collisions  of  ions  with  metal  surfaces  at  energies  between  1  au 
and  170  au.  A  set  of  localized  basis  functions,  consisting  of  generalized  Wannier 
functions  for  the  surface  and  s-  and  p-  atomic  functions  for  the  ion,  is  used  to 
separate  the  system  into  primary  and  secondary  regions.  An  effective  Hamiltonian 
and  time-dependent  equations  for  the  electron  density  matrix  are  obtained  in  the 
primary  region,  where  most  charge  transfer  occurs.  The  equations  for  the  electron 
density  matrix  are  solved  with  a  linearization  scheme.  The  method  is  suitable  to 


vii 


study  atomic  orbital  polarization  for  collisions  of  ions  and  surfaces.  A  model 
calculation  for  Na++W(110)  collisions  with  a  prescribed  trajectory  is  presented. 
The  interaction  potentials  between  the  W(110)  surface  and  Na  3  s  and  3p  orbitals 
are  calculated  from  a  Na  pseudopotential  and  a  step  potential  for  W(110).  Results 
show  that  the  yield  of  neutralized  atoms  oscillates  with  the  collision  energy  as 
the  result  of  the  near-resonance  charge  transfer  mechanism.  The  time-dependence 
of  the  density  matrix  provides  insight  on  the  dynamics  of  electron  transfer  along 
the  atomic  trajectory. 
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CHAPTER   1 
INTRODUCTION 


1.1  Collisional  Charge  Transfer  at  Surfaces 


Surfaces  and  various  physical  processes  occurring  near  surfaces  are  receiving 
increased  attention  in  present  scientific  research,  due  mainly  to  rapidly  growing 
applications  in  modern  analytical  techniques.  Researchers  in  fields  of  surface  sci- 
ence, chemical  physics,  solid  state  physics,  atomic  physics,  electrical  engineering, 
and  even  medical  science  are  carrying  on  extensive,  and  sometimes  overlapping, 
studies  and  analyses  involving  surfaces.  One  of  the  most  important  ways  to  un- 
derstand the  nature  of  surfaces  is  to  study  various  electronic  processes  occurring 
in  collisions  of  ions  with  them.  For  example,  one  can  learn  from  scattering  about 
the  topology  and  electronic  structures  of  surfaces,  amount  of  impurities  and  ad- 
sorbates  present,  surface  electron  and  lattice  relaxations,  and  interactions  between 
surfaces  and  atoms  and  molecules.  Particular  attention  has  been  paid  to  nona- 
diabatic  processes  accompanying  charge  transfer  between  scattered  particles  and 
surfaces  at  low  collision  energies  ranging  from  a  few  electron  volts  to  a  few 
thousand  electron  volts,  which  often  involve  strong  coupling  and  energy  transfer 
between  electronic,  translational,  vibrational,  and  rotational  degrees  of  freedom 
of  projectiles  and  surfaces. 

l 


Charge  transfer  in  ion-surface  collisions  may  lead  to  neutralization  or  ioniza- 
tion and  can  be  expressed  as 


E  +  A+  -*  E+  +  A 


and 


(1.1) 


E  +  A 


S±  +  A:F 


(1.2) 


respectively,  where  £  is  a  surface  such  as  W(110)  or  Ni(lll),  and  A  is  an  alkali 
atom  such  as  Li,  Na  or  K.  Similar  processes  occur  in  beam-foil  experiments  in 
which  charge  exchange  takes  place  when  ions  or  atoms  pass  through  thin  foils, 
instead  of  being  reflected  from  surfaces. 

The  two  main  mechanisms  through  which  charges  are  exchanged  are  the  near- 
resonance  process  and  the  Auger  process  [Hagstrum,  78].  In  a  near-resonance 
neutralization,  an  electron  in  the  conduction  band  of  the  solid  tunnels  to  an 
unoccupied  level  of  the  ion  which  lies  near  or  below  the  Fermi  level  of  the 
solid.  In  a  near-resonance  ionization,  an  electron  from  a  filled  excited  state  of  a 
neutral  atom  transits  to  an  empty  level  in  the  conduction  band  or  above  the  Fermi 
level  of  the  solid.  The  two  energy  levels  involved  in  these  processes  are  close, 
which  characterizes  the  near-resonance  processes.  The  Auger  process  involves 
two  or  more  electrons.  In  two-electron  Auger  neutralization,  one  electron  jumps 
from  the  metal  conduction  band  to  the  ion  ground  state  below  the  band  and,  to 
conserve  energy,  the  second  electron  is  excited  from  another  level  in  the  band  and 
may  leave  the  surface  and  be  detected  by  an  Auger  spectroscopy.  For  systems 


where  an  ion  affinity  level  is  in  resonance  with  the  conduction  band  of  the  metal, 
such  as  in  Na++W(110),  a  near-resonance  process  is  more  likely;  for  systems 
where  an  ion  ground  state  is  below  the  conduction  band,  such  as  H++Ni(lll),  the 
Auger  process  is  usually  the  dominant  charge  transfer  mechanism.  Complicated 
combinations  of  the  two  processes  are  also  possible  in  certain  circumstances. 

The  rise  of  interest  in  these  problems  has  promoted  an  extensive  effort  to 
develop  ultra-high  vacuum  (UHV)  and  modern  spectroscopy  techniques  to  observe 
and  detect  various  physical  properties  in  much  greater  detail  in  the  last  two  decades 
than  before  [Hagstrum,  77;  Hagstrum,  78].  In  these  spectroscopy  experiments, 
one  measures  the  intensity  of  scattering  products  (atoms,  photons  or  Auger 
electrons)  and  their  energy  and  angular  momentum  distributions  which  usually 
depend  on  the  energy  and  angle  of  incident  beams,  the  surface  properties  and, 
sometimes,  on  temperature  and  external  fields.  The  ion  scattering  spectroscopy 
(ISS)  [Smith,  71;  Hulpke,  75;  Taglauer  and  Hailand,  76;  Overbosch  et  al.,  80] 
measures  the  scattered  ion  intensity  as  a  function  of  incident  energy  and  reflected 
angle,  and  can  be  used  to  determine  the  topology  of  surfaces.  Oscillation  in  the 
intensity  of  neutralization  of  ions  in  the  low  energy  range  has  been  observed  for 
a  number  of  ion-surface  combinations,  and  is  ascribed  to  near-resonance  charge 
transfer  [Erickson  and  Smith,  75;  Rusch  and  Erickson,  77].  In  ion  neutralization 
spectroscopy  (INS)  experiments,  one  can  gain  knowledge  of  the  band  structure 
by  measuring  the  kinetic  energy  of  Auger  electrons  [Hagstrum,  75].  Time  of 
flight  (TOF)  experiments  detect  survival  probabilities  of  states  of  scattered  atoms 


and  their  energy  distributions  [Kasi  et  al.,  88].  So  far,  most  experiments  have 
been  done  on  noble  and  transition  metals,  simply  because  their  surfaces  are 
relatively  easy  to  clean  and  to  prepare.  However,  due  to  the  strong  interest 
in  electronics,  studies  on  semiconductor  and  other  material  surfaces  have  been 
carried  out  [Richard  and  Eschenbacher,  84]. 

In  atomic  collision  studies,  people  have  developed  new  instruments  which 
have  a  time  resolution  of  picoseconds  to  monitor  the  time  evolutions  of  inter- 
mediate atomic  states.  It  has  been  possible  to  measure  the  angular  distribution 
of  scattering  product  orbitals  [Hale  et  al.,  84;  Hertel  et  al,  85;  Andersen,  87; 
Campbell  and  Hertel,  87].  These  new  techniques  provide  valuable  information 
about  the  collision  dynamics,  and  it  is  expected  that  they  will  be  soon  used  in 
studies  of  ion-surface  collisions. 


1.2  A  Survey  of  Theoretical  Methods 

Advances  in  spectroscopy  techniques  in  the  last  two  decades  have  allowed  ex- 
tensive experimental  studies  of  ion-surface  collisions.  However,  the  development 
of  the  theory  lags  behind  that  of  experiment  in  this  field,  and  there  are  far  and 
away  more  experimental  data  than  can  be  understood.  This  is  because  the  systems 
we  are  dealing  with  are  very  challenging  ones  in  the  sense  that  they  lack  transla- 
tion symmetry.  Moreover,  they  involve  dynamical  many-electron  processes  with 
a  strong  and  time-dependent  coupling  between  the  electronic  degrees  of  freedom 
and  nuclear  motions.  The  traditional  methods  which  deal  with  systems  of  periodic 


lattices  or  of  a  few  electrons  are  either  not  suitable  or  have  to  be  modified.  It 
seems,  however,  that  theoretical  studies  are  accelerating  and  are  on  the  verge  of 
reproducing  some  of  the  experimental  data. 
1.2.1  Binary  Collision  Theory  of  Charge  Transfer 

In  this  theory  [Rusch  and  Erickson,  77;  Kasi  et  al.,  89]  the  scattering  of  an  ion 
from  a  surface  is  described  by  a  single  binary  collision  or  a  sequence  of  elastic 
collisions  between  the  ion  and  surface  atoms.  In  this  theory,  the  scattering  angle 
9  satisfies  the  relationship 


Ex 


cos9  +  (fi    —  sin  9) 


2M*/2 


(1.3) 


Eo  "  (1+fx)2^ 

where  mi  and  m2  are  the  masses  of  the  incidental  ion  and  surface  atom,  /z=m2/mi, 
and  Eo  and  E\  are  the  kinetic  energies  of  the  ion  before  and  after  the  scattering. 
The  yield  of  the  scattered  ions  for  a  given  scattering  angle,  Y(Eq,  9),  is 


Y(EQ,0)  oc  <j(Eo,0)  ■  P(E0,9) 


(1.4) 


where  <r(Eo,  0)  is  the  elastic  differential  cross  section  and  P(Eo,  9)  is  the  probability 
for  an  ion  to  remain  ionized  after  the  scattering.  In  classical  calculations  a  is  a 
monotonically  decreasing  function  of  the  incident  energy  and  P  is  given  by 


(1.5) 


P  oc  exp(—a/v±) 

where  a  is  a  constant  and  v±  is  the  ion  velocity  perpendicular  to  the  surface. 
This  theory,  in  which  the  yield  is  a  smooth  function  of  the  incident  ion  energy, 
fits  the  energy  dependence  of  4He+  scattered  from  Cu  and  certain  other  systems. 


However,  it  cannot  explain  the  oscillatory  charge  transfer  behavior  observed  in 
other  systems  which  is  believed  to  be  associated  with  nonadiabatic  charge  transfer 
mechanisms. 


1.2.2  Many-Electron  Treatments 

A  variety  of  theoretical  methods  [Tully,  77;  Brako  and  Newns,  81;  Holloway 
and  Gadzuk,  85;  Hood  et  al.,  85;  Lee  and  George,  85]  have  been  proposed  to  ac- 
count for  the  many-body  processes  accompanying  charge  exchange  in  ion-surface 
collisions.  The  most  widely  used  one  is  the  approach  based  on  the  Anderson- 
Newns  Hamiltonian  [Anderson,  61]  to  describe  the  ion-surface  scattering  [Blandin 
et  al.,  76;  Norskov  and  Lundqvist,  79;  Brako  and  Newns,  81;  Lang,  83].  In  this 
approach,  neglecting  spin  and  using  second  quantization  notations,  the  Hamilton- 
ian of  the  system  is  written  as 


.      H  =  eaC+Ca  +  ]T  ekC+Ck  +  J2  [VakCtCk  +  V:kC}Ca]  (M 

k  k 

where  C+,  Ca  are  the  creation  and  destruction  operators  corresponding  to  atomic 
state  xj>a,  C£,  Cjt  are  the  creation  and  destruction  operators  corresponding  to  solid 
state  0*,  ea  and  e*  are  the  energies  of  0a  and  0*  respectively,  and  Vak  =  {a\ V\k), 
where  V  is  the  perturbation  due  to  coupling  of  the  atom  to  the  metal.  The  last 
two  terms  describe  the  electron  hopping  between  tpa  and  ipk-  It  should  be  noted 
that  the  energy  of  the  atomic  level  is  time  dependent  through  the  ion's  trajectory 


RA(t),  that  is,  ea(t)  =  Haa  RA(t) 


A  set  of  equations  of  motion  of  electronic 


degrees  of  freedom  can  be  obtained  from  the  Heisenberg  equations  (with  h  =  1 ) 


.da 


i-^  =  -[H,  Ca]  =  ea(t)Ca  +  J2  KkWk 


(1.7) 


.dCk 
1  dt 


=  -[H,Ck]  =  ekCk  +  Vak(t)Ca 


(1.8) 


Eliminating  Ck,  it  yields  a  differential-integral  equation  for  Ca  which  can  be 
solved  by  numerical  procedures.  The  number  of  electrons  on  the  atom  after  the 
scattering  is  given  by 


MOO)  =  <V>a|C+(00)Ca(00)|V>a) 


(1.9) 


Most  studies  using  the  time-dependent  Anderson-Newns  Hamiltonian  have 
neglected  the  intra-atomic  Coulomb  repulsion  to  simplify  the  calculations.  For 
most  alkali-metal  collisions,  which  produce  few  negative  ions,  this  approximation 
seems  to  be  justified.  But  for  scatterings  like  H-W(llO),  there  is  a  significant 
fraction  of  H~  products  after  scattering,  and  the  intra-atomic  Coulomb  repulsion 
plays  an  important  role.  Instead  of  using  the  above  Hamiltonian  one  should  use 


ho 


+  E  [v*kC+Cka  +  V;kCtCaa]  +  U{t)na}nal 


(1.10) 


ha 


where  a  is  a  spin  index,  naT  =  C^Cah  nal  =  C+,Cai  ,  and  U(t)  is  the  effective 
intra-atomic  Coulomb  repulsion  which  depends  on  the  distance  between  the  ion 
and  the  surface.    Using  the  Hartree-Fock  approximation,  one  obtains  a  set  of 


differential-integral  equations  [Yoshimori  et  al.,  84;  Sulston  et  al.,  88]  which  can 
be  solved  numerically.  However,  the  qualitative  behavior  of  the  results  are  found 
to  be  not  much  different.  There  have  been  several  attempts  to  go  beyond  the 
Hartree-Fock  approximation  [Sebastian,  83;  85;  Kasai  and  Okiji,  87],  by  which 
some  preliminary  results  were  obtained. 


1.3  Existing  Problems 

Several  problems  inherent  in  the  theoretical  treatment  of  ion-surface  scattering 
will  be  addressed  here.  We  view  them  as  some  key  points  in  furthering  the 
theoretical  study  of  ion-surface  collisions. 

1.3.1  Electronic  States  and  Electronic  Couplings 

Ion-surface  systems  are  many  electron  systems  which  lack  translation  symme- 
try. This  brings  both  difficulty  and  challenge  to  theoretical  studies  of  the  problem. 
One  way  to  deal  with  surfaces  is  to  assume  that  the  electronic  wave  functions 
in  surfaces  are  the  same  as  those  of  an  infinite  periodic  system  and  to  treat  the 
colliding  ion  as  a  moving  defect  and  the  atomic  states  as  localized  states.  In 
this  way,  one  can  relatively  easily  calculate  surface  properties  and  the  scattering 
yield  [Brako  and  Newns,  81;  Lang,  83;  Shindo  and  Kawai,  86].  However,  this 
approximation  is  not  very  accurate,  and  usually  can  only  be  used  as  a  starting 
point.  At  surfaces,  as  the  results  of  broken  symmetry  and  the  strong  perturbation 
from  the  colliding  ion,  the  electronic  bands  are  distorted  and  the  electron  behavior 
is  different  from  that  of  bulks.  One  of  the  important  features  of  the  metal  and 


9 

semiconductor  electronic  structure  is  the  existence  of  the  localized  states.  While 

s-  and  p-electrons  are  generally  described  by  continuous  states,  d-electrons  are 
more  localized.  At  surfaces,  localized  states  can  also  be  created  by  impurities, 
adsorbates  or  chemisorbed  layers.  The  presence  of  localized  states  greatly  affects 
the  local  electronic  environment;  for  example,  an  absorbed  layer  of  alkali  on  metal 
surfaces  lowers  the  work  function  [Gomer,  75;  Medvedev  et  al.,  70]. 

Correctly  calculating  or  estimating  the  electronic  couplings  between  the  ion 
and  surface  is  obviously  important  in  theoretical  calculation.  But  current  theoreti- 
cal methods  do  not  provide  a  simple  and  yet  accurate  way  to  estimate  the  electron 
hopping  matrix  elements.  People  have  to  rely  on  semi-empirical  calculation.  For 
example,  in  the  time-dependent  Anderson-Newns  Hamiltonian,  Vak  measures  the 
coupling  between  the  atomic  states  and  surface  states  and  is  related  to  the  lifetime 
broadening  of  an  atomic  state,  A(r).  Assuming  that  k  and  t  dependence  of  V^ 
are  separated, 


V. 


ak 


RaW    =Vku(t) 


(1.11) 


where  V*  depends  only  on  k  and  «(t)  only  on  /,  one  has 


A(i)  =  A(C)M0I: 


(1.12) 


where  A(e)  is  defined  by 


(1.13) 


k 

which  can  be  approximately  evaluated  by  methods  developed  in  chemisorption 
theory.  But  a  common  practice  is  to  use  parameterization  of  A(f).  For  instance, 


10 

one  can  use  a  simple  exponential  form  which  is  independent  of  energy 


A(t)  =  A0exp(-7ZA)  (1.14) 

and  determine  parameters  Ao  and  7  from  experimental  data  [Brako  and  Newns, 
81;  Lang,  83;  Hood  et  al.,  85]  or  from  density  functional  calculation  [Lang  and 
Norskov,  83]. 

Because  of  the  electronic  coupling  from  the  surface  electrons,  the  atomic 
energy  ea  is  distance  dependent.  Again,  a  simple  function  form  is  often  used 
(with  the  electron  charge  e=l  au), 

£a(Za)  =  Ca(oo)  -  — -  (1.15) 


where  Z\  is  the  distance  between  the  ion  core  and  the  surface.  This  is  based  on 
the  approximation  that  the  change  of  ea  is  due  to  a  classical  image  charge. 

Although  the  parameterization  approach  is  widely  used,  it  ignores  the  fact 
that  the  electronic  coupling  is  strongly  dependent  on  the  rearrangement  of  charges 
during  collisions,  and  its  validity  at  short  distances  is  questionable.  Ideally,  the 
treatment  of  the  ion-surface  systems  should  contain  a  self-consistent  electronic 
structure  calculation  which  accounts  for  localized  states.  A  self-consistent  lo- 
calized orbital  method  has  been  developed  to  deal  with  transition  metal  surfaces 
and  chemisorption  on  metal  surfaces  [Smith  and  Gay,  75;  Smith  et  al.,  80].  The 
density  functional  method  [Hohenberg  and  Kohn,  64]  has  also  been  applied  to 
ion-surface  interactions  [Lang  and  Williams,  76]. 


1.3.2  The  Treatment  of  Extended  Systems 

From  the  theoretical  point  of  view,  it  is  easier  to  deal  with  either  a  system 
of  a  few  particles,  for  example,  small  molecules,  or  a  system  with  a  periodic 
structure,  such  as  single  crystals.  Ion-surface  systems,  on  the  other  hand,  are 
many-body  systems  which  lack  translation  symmetry.  To  combine  a  full  self- 
consistent  calculation  of  electronic  interaction  with  a  dynamic  description  of  ion- 
surface  collisions  is  a  difficult,  if  not  impractical,  task. 

In  developing  a  better  and  more  physical  description  of  ion-surface  systems 
a  very  important  fact  should  be  noticed,  namely  that  during  the  collisions,  ions 
interact  mainly  with  local  regions  of  the  surface,  while  the  remainder  of  the 
solid  is  relatively  unperturbed.     This  fact  provides  a  hint  that  it  is  possible 
to  properly  handle  electronic  motions  by  concentrating  on  these  local  regions. 
Olson  and  Garrison  [Olson  and  Garrison,  85]  have  used  a  cluster  in  place  of  a 
surface;  this  enables  them  to  employ  the  molecular  orbital  method  developed  in 
molecular  scattering.    Another  method  to  deal  with  surfaces  is  the  embedding 
technique  [Grimley  and  Mola,  74;  Kirtman  and  de  Melo,  81;  Feibelman,  85] 
used  in  chemisorption  studies,  which  uses  a  higher  level  treatment  of  electronic 
interactions  within  a  molecular  complex  (a  small  region  of  the  surface  plus 
adsorbates)  and  a  simple  description  outside  of  the  complex,  de  Melo  et  al.  [de 
Melo  et  al.,  87]  have  developed  a  self-consistent  method  using  a  density  matrix 
and  applied  it  to  the  Anderson-Newns  Hamiltonian  for  a  one  dimensional  system. 
McDowell  [McDowell,  85]  uses  an  embedding  technique  to  obtain  generalized 
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Langevin  type  equations  for  the  spin  orbitals  in  a  primary  zone  which  couples 
with  the  secondary  zone  through  driving  terms.  These  treatments  show  that, 
by  concentrating  on  a  small  number  of  orbitals,  it  is  possible  to  achieve  self- 
consistency.  Another  common  feature  in  these  treatments  is  the  use  of  localized 
spin  orbitals,  which  allows  one  to  naturally  deal  with  localized  phenomena  [Feng 


et  al.,  91]. 

1.3.3  Effect  of  Nuclear  Motion 


While  the  electronic  motions  are  treated  quantum  mechanically,  it  is  reason- 
able to  assume  that  the  nuclear  motion  evolves  according  to  the  classical  mechan- 
ics for  the  collision  energy  in  the  hyperthermal  range.  In  this  classical-quantal 
approach  the  nuclear  motion,  RA(t),  can  be  treated  at  several  levels.  In  the  simple 
classical  treatment,  the  nuclear  trajectory  is  fully  determined  by  a  classical  ion- 
surface  potential.  In  an  improved  semiclassical  treatment,  the  nuclear  potential  is 
coupled  with  the  electronic  motion.  The  classical  trajectory  is  valid  only  in  the 
high  energy  range  [Tully,  76]  as  demonstrated  in  gas-phase  collisions.  The  reason 
is  that  at  lower  energies,  trajectories  are  sensitive  to  the  detail  of  chemical  inter- 
actions among  the  atoms  and  depend  on  the  electronic  states.  It  has  been  found 
that  in  H-p  collisions,  the  classical  trajectory  would  give  rise  to  a  significant  error 
when  the  collisional  energy  is  lower  than  100  ev  [Runge  et  al.,  90]. 

Looking  from  another  angle,  the  nuclear  motion  can  be  treated  either  by 
the  so-called  trajectory  approximation  or  by  a  multichannel  procedure.  In  the 
trajectory  treatment,  the  nuclear  trajectory  is  uniquely  defined  by  a  single  fixed 
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potential,  which  can  either  be  a  pure  classical  one  or  an  effective  potential  which 
contains  the  coupling  with  the  electronic  degrees  of  motion.  One  example  is  the 
straight  line  trajectory  resulted  from  a  hard  wall  potential.  Runge  et  al.  [Runge 
et  al.,  90]  determine  the  ion  trajectory  by  an  effective  force  which  is  dependent 
on  the  gradient  of  the  electronic  density  matrix.  It  has  instead  been  suggested 
using  a  trajectory  approximation,  constructing  multidimensional  potential  energy 
hypersurfaces  describing  various  atomic  configurations  and  allowing  trajectories 
to  hop  back  and  forth  between  hypersurfaces  [Tully,  77;  Holloway  and  Gadzuk, 
85;  Newns,  85]. 


A  useful  procedure  in  dealing  with  molecule  collisions  is  the  eikonal  method 
developed  by  Micha  [Micha,  83].  In  this  method  the  nuclear  variables  are  coupled 
with  the  electronic  ones  and  both  must  be  found  self -consistently.  An  application 
of  the  eikonal  method  to  ion-surface  collision  is  by  Olson  and  Garrison  [Olson 
and  Garrison,  85]  who  model  the  surface  by  a  small  cluster. 

1.3.4  Atomic  Orbital  Polarization 

In  atomic  scattering,  the  scattered  atoms  can  be  in  different  electronic  states 
and  their  distribution  of  electronic  angular  momentum  can  be  anisotropic.  This 
causes  orbital  alignment  and  orientation  [Hippler,  85].  Experimentally,  this 
requires  a  partial  or  full  determination  of  the  atomic  states  after  the  collision  in 
contrast  to  the  conventional  study  which  measures  the  differential  cross  section. 
Orbital  polarization  has  been  the  subject  of  much  theoretical  work  [Fano  and 
Macek,  73;  Andersen  and  Nielsen,  87;  Nielsen  and  Andersen,  87],  in  which  the 
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density  matrix  [Blum,  81]  is  extensively  used.  Although  the  use  of  the  density 
matrix  in  scattering  theory  is  not  new,  the  study  of  orbital  polarization  provides 
a  much  more  severe  test  of  a  theory  than  does  the  study  of  cross  section,  this 
is  because  not  only  the  diagonal  elements,  but  also  the  off-diagonal  elements 
of  the  density  matrix  are  required  to  determine  the  alignment  and  orientation 
parameters.  Of  course,  the  study  of  orbital  polarization  will  lead  to  a  deeper 
insight  into  the  interaction  and  mechanisms  involved  in  the  collisions.  Until 
now,  no  detailed  work  on  orbital  polarization  in  ion-surface  collision  had  been 
reported  either  experimentally  or  theoretically.  It  is  just  a  matter  of  time  before 
the  experimental  techniques  and  theory  are  developed  in  this  field. 

1.4  Outline  of  the  Chapters 

This  dissertation  presents  a  theoretical  study  of  charge  transfer  in  ion-surface 
collisions  based  on  the  time-dependent  molecular  orbital  approach.  A  partition- 
ing technique  is  used  to  divide  a  ion-surface  system  into  a  primary  region  which 
contains  a  few  centers  on  the  surface  and  the  scattering  ion,  and  a  secondary 
region  containing  the  remainder  of  the  surface.  The  charge  transfer  is  calculated 
by  solving  the  time-dependent  Hartree-Fock  (TDHF)  equation  for  the  electron 
density  matrix  in  the  primary  region  which  couples  with  the  secondary  region 
electronically.  This  approach  permits  determination  of  the  atomic  states  during 
and  after  the  scattering  and  calculation  of  the  alignment  and  orientation  parame- 
ters. In  the  following  chapters  atomic  units  are  used;  therefore,  electron  charge 
e=l,  electron  mass  me=\  and  h=\. 
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In  Chapter  2,  we  present  the  basic  framework  of  our  approach  to  the  ion- 
surface  collisions.  TDHF  equations  are  used  to  describe  the  time  evolution  of  the 
molecular  orbitals.  The  formalism  is  constructed  in  terms  of  an  electronic  density 
matrix  because  it  can  simplify  computation  and  is  convenient  in  the  calculation  of 
alignment  and  orientation  parameters.  In  the  later  part  of  the  chapter,  the  TDHF 
equation  is  rewritten  in  the  basis  of  the  travelling  atomic  orbitals  (TOA). 


In  Chapter  3,  we  define  the  parameters  for  average  electronic  population, 
electronic  multipoles  and  polarization  of  atomic  orbitals  by  introducing  tensor 
operators  and  irreducible  operators.  The  orientation  and  alignment  parameters 
can  be  calculated  from  the  density  matrix.  The  orbital  polarization  parameters  for 
a  subsystem  are  also  expressed  by  the  density  matrix. 

In  Chapter  4,  we  present  a  local  time  linearization  procedure  for  solving 
TDHF  equations.  The  method  is  based  on  the  assumption  that,  in  a  very  short 
time  interval,  the  solution  of  the  TDHF  equation  is  a  linear  perturbation  caused 
by  the  motion  of  the  atomic  core,  added  to  correct  the  evolution  of  the  system 
first  calculated  as  if  the  nuclei  were  fixed.  Using  this  procedure,  the  linearized 
equations  for  the  perturbations,  which  contain  a  time-dependent  driving  term 
due  to  the  nuclear  motions,  are  constructed  and  solved  in  conjunction  with  the 
equations  for  fixed  nuclear  positions.  The  procedure  is  presented  for  a  general 
case.  In  the  later  part  of  the  Chapter,  we  apply  this  procedure  to  a  simple  case  in 
which  the  electron-electron  interaction  is  neglected  and  obtain  analytical  solutions 
for  the  linearized  TDHF  equations. 
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In  Chapter  5,  a  partitioning  procedure  for  ion-surface  collisions  is  presented. 

The  system  is  divided  into  a  primary  region  which  contains  the  ion  and  a  few 
centers  on  the  surface  and  a  secondary  region  which  is  the  remainder  of  the 
surface.  The  density  matrix  is  approximated  in  the  secondary  region  by  an 
unperturbed  one,  which  is  justified  when  the  collisional  energy  of  the  ion  is 
not  very  low.  The  effective  equations  in  the  primary  region  are  coupled  with 
the  secondary  region  through  a  driving  term.  To  fit  the  asymptotic  behavior  of 
the  system  after  the  partitioning,  a  correction  term  is  added  to  the  partitioned 
Hamiltonian.  Combining  the  linearization  procedure  described  in  Chapter  4,  we 
obtain  the  effective  equations  in  the  primary  region  and  their  formal  solutions. 
Again,  we  further  illustrate  this  partition  method  by  applying  it  to  a  special  case 
when  the  electron  interaction  is  neglected. 

In  Chapter  6,  we  construct  a  set  of  localized  basis  functions  to  be  used  in  the 
time-dependent  molecular  orbital  calculation.  For  the  surface  part,  we  introduce  a 
set  of  generalized  Wannier  functions  (GWFs)  which  are  orthonormal  and  localized 
at  centers  of  the  surface.  A  variational  procedure  is  employed  to  generate  these 
functions.  As  an  example,  we  calculate  the  generalized  Wannier  functions  for 
a  three  dimensional  jellium  solid.  The  basis  functions  associated  with  the  ion 
core  are  assumed  to  be  atomic  orbitals  and  are  obtained  from  a  pseudopotential 
for  the  ion  core.  Our  basis  functions  are  time-dependent  and  in  the  form  of 
linear  combinations  of  Gaussians.  The  overlap  matrix  and  Hamiltonian  matrix 
are  constructed  in  this  basis  set. 
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In  Chapter  7,  we  discuss  the  computational  aspects.  At  first,  we  present  our 
algorithm  for  numerical  integration  of  the  time-dependent  Hartree-Fock  equations 
as  well  as  the  flowchart  of  the  computational  program.  Then  the  stability  and 
the  convergence  of  the  numerical  integration  are  discussed  at  length.  Effects  of 
some  computational  parameters  are  investigated,  including  step  size,  initial  and 
final  positions  of  the  ion,  and  tolerances,  among  others. 

In  Chapter  8,  we  present  an  application  of  this  approach  and  its  results  to 
the  neutralization  in  a  normal  collision  between  a  Sodium  ion  and  a  Tungsten 
(110)  surface.  An  atom-surface  interaction  potential  is  constructed  to  determine  a 
prescribed  trajectory.  By  using  methods  described  in  previous  chapters,  we  obtain 
the  time  evolution  of  the  density  matrix  from  which  the  electronic  population  and 
orbital  polarization  parameters  are  calculated.  We  also  study  the  effects  of  the 
collision  energy  on  the  final  yield  of  the  neutral  atoms  and  their  polarization  after 
the  collision,  and  compare  the  results  with  experimental  data.  We  also  briefly 
discuss  the  application  of  the  approach  to  other  systems,  and  to  collisions  at  other 
scattering  angles. 

In  Chapter  9,  we  discuss  the  features  of  our  time-dependent  molecular  method 
and  its  applications  to  charge  transfer  in  atom-surface  collisions  in  hyperthermal 
energy  range.  Analysis  and  conclusions  are  related  to  the  physical  features  and 
comparison  with  other  theoretical  methods  in  the  field.  We  also  discuss  the 
approximations  used  in  the  present  study.  Finally,  we  offer  suggestions  for  future 
research. 
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In  Appendix  A,  we  calculate  matrix  Bx  which  is  the  expansion  coefficient 

matrix  of  the  Wannier  functions  for  free  electrons  in  terms  of  Gaussians. 

In  Appendix  B,  we  calculate  the  total  band  energy  of  a  jellium  slab  using  the 

generalized  Wannier  functions. 

In  Appendix  C,  we  show  how  the  Fermi  energy  is  calculated  for  a  finite 
jellium  slab. 


CHAPTER  2 
THE  TIME-DEPENDENT  HARTREE-FOCK  (TDHF)  APPROACH 


The  time-dependent  Hartree-Fock  approach  is  suitable  for  studying  the  many- 
body  collision  dynamics  of  our  problem.  Using  a  one-electron  effective  field  it 
can  describe  the  evolution  of  collisional  systems  and  calculate  dynamic  parameters 
to  compare  with  experimental  data.  One  of  the  advantages  of  the  TDHF  method 
is  that  the  self-consistency  between  the  effective  field  and  the  orbital  coefficients 
or  density  matrix  is  achieved  automatically  by  solving  time-dependent  differen- 
tial equations  without  having  to  perform  the  self-consistent  iteration  procedure 
required  in  the  time-independent  case.  TDHF  has  been  extensively  applied  to 
collision  problems  in  nuclear  physics  [Ring  and  Schuck,  80;  Negele,  82;  Davies 
et  al.,  82],  and  recently  applications  in  atomic  molecular  collisions  have  been 
reported  [Kulander  et  al.,  82;  Devi  and  Garcia,  83;  Gazdy  and  Micha,  86;  87]. 
Micha  and  Gazdy  [Micha  and  Gazdy,  87]  have  proposed  a  variational  procedure 
which  improves  the  accuracy  of  transition  amplitude  calculations  using  TDHF 
trial  functions.  Yoshimori  et  al.  [Yoshimori  et  al.,  84]  have  discussed  using  the 
TDHF  method  in  charge  transfer  in  ion-surface  collisions. 

In  this  Chapter,  we  apply  the  TDHF  method  to  charge  transfer  in  ion-surface 
collisions.     In  Sect.     2.1,  by  introducing  time-dependent  molecular  orbitals, 
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we  derive  the  TDHF  equations  for  the  electron  density  matrix.  The  reason 
we  use  the  density  matrix  [McWeeny,  60;  Cohen  and  Frishberg,  76]  is  that 
it  simplifies  calculation  for  a  large  system.  Another  advantage  is  that  many 
interested  properties  of  collisions  can  be  defined  and  readily  calculated  from  it. 
One  of  the  key  problems  in  applying  the  TDHF  approach  to  ion-surface  studies 
is  choosing  the  basis  set.  With  the  great  number  of  electrons  in  the  systems, 
choosing  a  complete  basis  set  would  make  any  practical  calculation  out  of  the 
question.  One  has  to  truncate  the  complete  basis  set  in  a  proper  way  so  that  the 
introduced  error  is  minimal.  However,  we  will  leave  this  problem  to  Chapter  5. 
In  this  Chapter,  we  first  describe  the  ion-surface  system.  In  Sect.  2.2  we  derive 
the  basic  formalism  without  specifying  the  particular  form  of  the  basis  set  except 
for  the  condition  that  each  basis  function  is  associated  with  a  certain  center  of 
the  system.  In  Sect.  2.3  the  TDHF  equation  is  rewritten  in  the  basis  of  travelling 
atomic  orbitals  (TAO). 

2.1  The  Atom-Surface  System 

Let's  consider  a  system  of  a  scattering  ion  and  a  semi-infinite  solid.  We  will 
assume  the  solid  has  a  conduction  band  of  continuous  one-electron  states  and 
some  localized  states  (This  is  in  contrast  to  some  approaches  which  are  limited 
to  continuous  states  only).  For  the  colliding  atom  we  assume  it  has  a  core  and 
several  valence  levels.  Because  the  deeper  levels  of  the  solid  and  the  core  states 
of  the  ion  are  tightly  bound  and  usually  do  not  participate  in  charge  transfer  in 
the  energy  range  we  are  interested  in,  they  are  not  considered.    For  the  same 
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reason,  highly  excited  atomic  states  and  those  solid  states  highly  above  the  Fermi 
level  are  not  considered  either.  We  also  assume  that  the  atomic  levels  are  in  the 
range  of  or  close  to  the  conduction  band  and  that  the  near-resonance  tunnelling 
of  electrons  is  the  dominant  charge  transfer  process.  The  thermal  motion  of  the 
nuclei  in  the  solid  can  be  neglected  at  the  collision  energies  of  interest,  so  that 
the  only  moving  nucleus  is  that  of  the  scattering  ion.  For  simplicity,  we  assume 
that  the  collision  energy  is  high  enough  that  the  trajectory  approximation  can  be 
used.  In  this  approximation,  the  motion  of  the  ion  nucleus  can  be  described  by 
a  known  function  Rj^t)  determined  by  an  effective  interaction  potential  between 
the  ion  core  and  the  surface.   In  Chapter  9  we  will  discuss  the  validity  of  this 
approximation  and  the  possible  improvements  which  are  necessary  at  low  collision 
energies. 


2.2  TDHF  Equations  of  Density  Matrix 

Let  ^(t)  be  the  electronic  wavefunction  which  satisfies  the  time-dependent 
Shrodinger  equation 


d 


(2.1) 


where  H(t)  is  the  total  Hamiltonian  of  the  system.    Using  the  time-dependent 
Hartree-Fock  approximation,  the  TDHF  wavefucntions  have  the  form  of 


v(t)  =  Xd'teM*j,*)] 


(2.2) 
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where  Nei  is  the  number  of  electrons  and  ipi(x,  t)  is  the  ith  spin  orbital  for  electron 
variables  x  =  (r, (,).   We  write  tpi(x,t)  as 


&■(*,*)  =  tf(r,<)7.-(0 


(2.3) 


where  <j>J(r,t)  is  the  /th  time-dependent  molecular  orbitals  for  spin  7  and  7i(() 
is  the  corresponding  spin  function.  The  electron  density  operator  p  is  given  by 


(2.4) 
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where  the  summation  is  over  the  occupied  orbitals,  which  is  assumed  to  satisfy 
the  TDHF  equation  [Dirac,  30] 


ih%.pt  =  Ftp1  -  p*f* 
ot 


(2.5) 


where  F  is  the  Fock  operator  which  can  be  written  as  [Pople  and  Beveridge,  70] 


where 


rr{xuiy=£(?ut)+&(xi,t) 


Wut)  m  -1  V?  +VA(ri,t)  +  VM{rut) 


(2.6) 


(2.7) 


is  the  core  Hamiltonian  operator  with  Va  the  potential  from  the  atomic  core  and 
Vm  the  potential  from  atomic  cores  in  the  solid;  and 

fr(xut)  =  J2Jdx2^(x2,t)^-[l  -Vl2]ipi(x2,t)  (2.8) 

is  the  electron  Coulomb  potential  operator,  where  V\2  is  the  permutation  operator 
exchanging  electrons  only  between  spin-orbitals  of  spin  7.  Introducing  a  set  of 
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time-dependent  basis  functions  {^(r,i)},  fx=l,  2,  .  .  .  Nb,  the  molecular  orbital 
is  written  as  a  linear  combination  of  the  basis  functions  so  that 

tf  (r,  t)  =  J2  cl(t)M?>  *)  »  =  1,  2,  •  •  •  (2.9) 

ft 

where  c] At)  is  a  time-dependent  coefficient.  The  density  matrix  P  is  defined  as 

f  =  I0P7(£I  (2-10) 

where  |£)  and  (£|  are  row  and  column  matrices  of  £«  orbitals,  respectively.  The 
matrix  element  of  P  is  [Pople  and  Beveridge,  70;  Szabo  and  Ostlund,  82] 


i£occ 


The  Fock  matrix  F^  is  defined  as 

Inserting  equation  (2.6)  into  the  above  definition  we  have 

F7  =  H  +  G7(p7,P7') 


(2.11) 


(2.12) 


where 


H  =  K  +  VA  +  VM 


(2.13) 


(2.14) 


is  the  core  Hamiltonian  matrix,  G^  the  Hartree-Fock  electron-electron  interaction 
matrix,  Va  the  atomic  potential  matrix  and  Vm  the  surface  potential  matrix.  The 
matrix  elements  of  F^  are 
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where  7'  is  the  spin  opposite  to  7  and 


(fiu\\\a)  =  ((ii/\\<t)  —  (fxa\\u) 


(2.16) 


with 


(HAcr)  =  / 


73 /3„_  £*(-? 


1 


n2 


(2.17) 


To  derive  the  TDHF  equation  for  the  density  matrix,  let's  define  a  matrix 

d 


n  = 


®vfl 


(2.18) 


Substitute  (2.10)  into  (2.5),  the  left  hand  side  is 


and  the  right  hand  side  is 


IO)P7(^l  +  lo|pT(ei  +  IOP7 


fM 


(2.19) 


fv  -  pfr  =  nop^ti  -  iop7<^7 


(2.20) 


Multiply  this  equation  by  (£|  from  the  left  and  by  |£)  from  the  right,  we  have 


i«P7S  +  iSP7S  +  iSP7&  =  F7P7S  -  SP7F7 


(2.21) 


where  S  =  (<f|<f)  is  the  overlap  matrix.  Multiply  both  sides  of  the  above  equation 
by  S_1  we  obtain  the  time-dependent  Hartree-Fock  equation  for  the  density  matrix 


jP7  =  S-^P7  -  P^S"1  -  iS-^P7  -  zP^S-1 


(2.22) 


The  last  two  terms  on  the  right  hand  side  of  Eq.    (2.22)  arise  from  the  time- 
dependence  of  the  basis  functions. 
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2.3  TDMOs  as  Linear  Combinations  of  Travelling  Atomic  Orbitals 


When  the  molecular  orbitals  are  expanded  in  a  basis  set  of  atomic  functions, 
it  sometimes  introduces  artificial  couplings  at  large  distance  which  originates  in 
the  dependence  of  the  atomic  orbitals  on  the  position  of  the  moving  nuclei  [Bates 
and  McCarroll,  1958].  To  avoid  this  one  can  choose  the  atomic  basis  functions 
in  the  form  of  travelling  atomic  orbitals  (TAO).  We  choose  £«  to  be  a  travelling 
atomic  orbital  (TAO)  associated  with  the  rath  center  at  Rm(t), 


tii(r,t)  =  Xn  f-Rm(t)  Tm(r,t) 


(2.23) 


where  f  is  the  position  vector  with  respect  to  the  origin  of  the  reference  frame, 
Xn  is  the  atomic  orbital  centered  at  Rm(t)  and  the  translation  factor  Tm(f,  t)  is 


defined  by 


Tm{f,t)  =  exp 


imr 


(2.24) 


where  me  is  the  electron  mass  and  vm  =  dRm/dt  is  the  local  velocity  of  the  rath 
center  in  a  space-fixed  system.   Thus, 


R'W-CSjj6***-** 


) 


and 


iW«(&|gj6,) 

d 


) 


(6.l(gj)„W  +  (&l(«.-v)w 


Using  (2.22)  we  find 


i-   |2 


Qpu  =  —im( 


(2.25) 


(2.26) 


l^i-5^  +  (Xtl\T*mTn(vn  ■  Vn)|x„)  (2.27) 
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where 


(X»\KTn(vn.Vn)\Xl,) 

=  J  d3rx;(f)K(r)Tn(r)[[vn  •  Vn)x,(0 

— *  — » 

Noticing  that  V|x«/)  =  -V„|xi/),  the  kinetic  energy  matrix  in  this  basis  is 

=  ^jpU,,, + i<xiiister.(*  •  vB)ix„>  +  (xniA'M   (2-29) 

where  Eq.    (2.26)  has  been  used  and 


Substituting  Eq.   (2.30)  into  Eq.   (2.22)  it  becomes 


with  a  modified  Fock-like  matrix 


iP7  =  S_1F7P7  -  P7F7tS_1 


K„  =  (Xll\K\Xu)  =  jd\  T^(r,t)Tn(r,t)  x^t)KX^t)  (2. 


30) 


(2.31) 


F7  =  H  +  G" 


where 


H  =  K  +  Vi4  +  VM 


(2.32) 


(2.33) 


Equation  (2.31)  is  our  basic  equation.  It  appears  in  a  simpler  form  in  the 
TAO  basis,  but  the  price  paid  is  that  the  matrix  F7  is  no  longer  Hermitian. 

This  equation  can  be  solved  numerically  if  the  nuclear  trajectory  is  known, 
which  may  or  may  not  include  the  effect  of  the  coupling  between  the  electron 
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and  nuclear  degrees  of  freedom.  In  the  eikonal  treatment,  this  equation  should 
be  solved  with  the  equations  of  motion  for  nuclei  simultaneously.  Runge  et  al. 
[Runge  et  al.,  90]  have  applied  this  method  to  charge  transfer  in  the  collision 
between  a  hydrogen  atom  and  a  proton. 

So  far,  we  have  put  no  restrictions  on  the  choice  of  basis  functions  except 
that  the  functions  are  localized  at  certain  centers  of  the  system.  In  Chapter  6  we 
will  discuss  in  detail  the  choice  of  basis  functions,  which  is  a  key  part  to  our 
approach  for  handling  extended  systems. 


' 


CHAPTER  3 

AVERAGE  ELECTRONIC  POPULATIONS,  ELECTRIC 

MULTIPOLES  AND  ORBITAL  POLARIZATION 


One  of  the  parameters  which  characterizes  scattering  phenomena  is  the  cross 
section,  which  has  the  dimension  of  area  and,  roughly  speaking,  measures  the  size 
of  electron  clouds.  The  total  cross  section  and  the  differential  cross  section  for 
state  to  state  transitions  can  be  measured  in  experiments  and  provide  fundamental 
information  about  excitation  mechanisms.  Thus  extensive  efforts  have  been  made 
to  develop  theories  and  methods  to  calculate  the  cross  section.  On  the  other  hand, 
during  collisions  with  targets,  absorption  of  energy  will  excite  atoms  to  various 
states  which  can  be  anisotropic.  Electron  clouds  not  only  change  their  sizes  but 
also  change  their  shapes  and  rotations,  which  is  termed  orbital  polarization.  The 
shape  change  and  rotation  of  the  electron  cloud  are  characterized  by  alignment  and 
orientation.   Recent  advances  in  experimental  techniques  have  made  it  possible 
to  completely  determine  the  states  of  scattered  atoms.  For  example,  in  the  H+-H 
experiment  [Hippler  et  al.,  86],  H  experiences  a  2s  to  2p  excitation  and  the 
distribution  of  electrons  on  m=\  and  m=-\  states  varies  with  the  collision  energy. 
The  orbital  polarization  reflects  some  of  the  specific  and  subtle  aspects  of  the 
collision  processes  which  can  not  be  fully  unveiled  in  the  cross  section  study. 
Obviously,  knowledge  of  orbital  polarization  properties  can  offer  a  deeper  insight 
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into  the  collision  mechanism  and  a  more  detailed  understanding  of  the  dynamics 
of  collisions.    Such  a  study  will  also  provide  a  sensitive  test  to  the  scattering 
theories  and  models. 

The  density  matrix  method  has  been  proved  to  be  a  powerful  tool  in  inves- 
tigations of  orbital  polarization  in  scatterings  [Fano  and  Macek,  73;  Andersen, 
87;  Hippler  et  al.,  86].   Using  the  language  of  the  density  matrix,  the  diagonal 
elements  are  related  to  the  cross  section,  and  the  off-diagonal  elements  contain 
the  information  about  the  shape  and  rotation  of  the  electron  cloud.  Usually,  the 
density  matrix  is  constructed  to  calculate  the  orbital  polarization  parameters  for 
the  electron  cloud  associated  with  the  scattered  atom.  In  this  way  one  assumes 
that  the  electron  cloud  contains  only  the  electrons  in  the  orbitals  of  the  scattered 
atoms  and  that  the  contribution  from  the  target  electrons  is  excluded.   Such  an 
approach  is  suitable  for  most  experimental  situations  where  the  scattered  atoms 
are  detected  far  away  from  the  target  but  before  or  in  the  course  of  decaying 
from  their  excited  states.    However,  at  short  or  medium  distances,  the  atomic 
orbitals  are  mixed  with  the  surface  orbitals  and  the  electron  cloud  associated  with 
the  atom  contains  also  the  contribution  from  the  surface  electrons.   The  above 
approach  does  not  seem  satisfactory  from  our  theoretical  point  of  view,  since  we 
are  interested  in  following  the  time  evolution  of  the  electron  cloud  which  requires 
to  describe  the  orbital  polarization  of  the  atom  at  all  distances.  The  definitions  of 
the  orbital  polarization  parameters  should  be  consistent  with  that  for  an  isolated 
atom  as  the  distance  between  the  atom  and  the  surface  becomes  infinitely  large. 
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In  this  Chapter,  we  use  the  electronic  density  matrix  to  define  parameters 
characterizing  the  orbital  polarization.  Unlike  other  studies  [Fano  and  Macek,  73; 
Nielsen  and  Andersen,  87],  we  start  with  the  electronic  density  matrix  of  the  full 
system,  i.e.,  the  atom  plus  the  surface,  and  then  define  the  orbital  parameters  for 
a  subsystem  which  either  is  the  scattered  atom  alone  or  is  a  molecular  complex 
containing  the  atom  and  a  part  of  the  surface.  Because  the  mixing  between  the 
atomic  orbitals  and  the  surface  orbitals  are  taken  into  account,  the  electron  cloud 
associated  with  the  scattered  atom  contains  also  the  contribution  from  the  surface 
electrons.  This  permits  us  to  examine  the  evolution  of  the  electron  cloud  and 
gives  a  dynamic  picture  of  the  orbital  polarization  at  all  the  distances  during 
the  collision.  As  the  distance  between  the  atom  and  the  surface  becomes  very 
large,  the  contribution  from  the  surface  electrons  becomes  insignificant  and  our 
polarization  parameters  for  the  atom  orbitals  are  asymptotically  equal  to  those 
defined  in  other  approaches. 

In  Section  3.1,  we  describe  two  coordinate  systems:  the  collision  frame  and 
the  natural  frame.  The  collision  frame  seems  suitable  to  describe  the  scattering;  the 
natural  frame  is  more  convenient  to  picture  orbital  polarization.  In  Section  3.2, 
using  coordinate  tensor  operators,  we  define  electric  multipoles  which  provide 
a  picture  of  the  spacial  distribution  of  the  electron  cloud.  In  Section  3.3,  the 
alignment  and  orientation  parameters  are  defined  by  use  of  the  irreducible  tensor 
operators  and  expressed  them  in  terms  of  the  density  matrix.  In  Section  3.4, 
the  polarization  parameters  are  defined  for  a  subsystem  which  contains  either  the 
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scattered  atom  alone  or  the  atom  plus  a  small  part  of  the  surface.  This  is  useful  for 
comparison  with  experimental  data  from  spectroscopies.  The  connection  between 
the  orbital  polarization  parameters  for  a  subsystem  and  that  for  an  isolated  atom 
is  discussed. 

3.1  Coordinate  Frames 



Since  symmetry  plays  an  essential  role  in  the  study  of  scatterings  and  orbital 
polarization,  it  is  the  first  thing  one  should  look  at  when  choosing  a  coordinate 
frame.  Let's  consider  the  symmetry  of  the  p-orbitals.  The  three  eigenfunctions 
of  the  angular  momentum  operator  for  7=1,  |p+i),  |p0)  and  |p_x)  correspond  to 
the  magnetic  quantum  number  m=\,  0  and  -1  respectively.  We  can  also  use  three 
real  p-orbital  functions  |px),  |py)  and  |pz)  which  are  symmetric  along  x-,  y-  and 
z-axes  respectively.  The  orbital  |pz)  =  |p0)  has  a  negative  reflection  symmetry 
with  respect  to  the  X-Y  plan  and  others  have  positive  reflection  symmetry. 

The  collision  plane  is  determined  by  the  incoming  velocity  vector  vm  and 
outgoing  velocity  vector  vout  of  the  projectile.  For  an  atom-surface  collision,  the 
collision  plane  is  perpendicular  to  the  solid  surface.  The  collision  frame  (Xc, 
Yc,  Zc)  is  defined  such  that  the  Xc-Zc  plane  coincides  with  the  collision  plane 
with  the  Zc  axis  perpendicular  to  the  surface  and  that  the  Xc  axis  parallel  to  the 
surface,  and  the  Xc-Yc  plane  coincides  with  the  solid  surface  with  the  Yc  axis 
perpendicular  to  the  collision  plane,  see  Fig.  1.  Because  many  collision  systems 
have  certain  type  of  symmetry  about  the  Zc  axis  or  with  respect  to  the  collision 
plane,  this  frame  is  convenient  for  collision  problems  [Andersen,  86]. 


*c  '  YN 


Figure  1.1.  The  coordinate  frames  used  for  the  description  of  orbital 
polarization  of  the  electron  cloud.  The  collision  frame  (Xc,  Yc,  ZJ  and 
the  natural  frame  (Xn,  Yn,  Zn)  are  shown.  The  incoming  velocity  vector 
v[n  and  the  outgoing  velocity  vector  vout  are  in  the  collision  plane  which 
is  perpendicular  to  the  solid  surface  and  coincides  with  the  Xc-Zc  plane  of 
the  collision  frame  and  the  Xn-Yn  plane  of  the  natural  frame. 
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In  the  natural  frame  (X„,  Yn,  Zn),  as  shown  in  Fig.  1.1,  the  Xn-Yn  plane 
is  in  the  collision  plane  with  the  Xn  axis  perpendicular  to  the  solid  surface,  the 
Yn-Zn  plane  coincides  with  the  solid  surface  with  the  Z„  axis  perpendicular  to  the 
collision  plane.  This  frame  seems  to  be  "natural"  to  the  analysis  of  the  angular 
momentum  transfer  and  orbital  polarization,  since  the  expectation  value  of  the 
z-component  of  angular  momentum  is  related  to  the  orientation  and  alignment 
[Andersen,  86]. 

3.2  Electric  Multipoles 

One  way  to  describe  the  shape  of  the  electron  cloud  is  to  consider  its  electric 
multipoles.  For  this  purpose  we  introduce  a  set  of  tensor  operators  in  the  Cartesian 

coordinates, 

7<°>  =  1 


/(?)  =r2c.._o 


(3.1) 


where  i,  j=x,  y,  z,  and  the  super  indices  in  the  parentheses  indicate  the  ranks  of 
the  tensor  operators.  The  electric  multipoles  are  defined  as 


P(k)  _  p(k)) 


(3.2) 


where  I(k)  represents  the  tensor  operator  defined  in  Equation  (3.1),  and  the  symbol 
<  >  indicates  the  ensemble  average  which  can  be  expressed  in  terms  of  the  density 
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(3.3) 


operator  p  and  the  electronic  density  matrix  PT  in  a  basis  set  {^} 

/T(fc)\  =       V  /    =    1    "v  'V]i 

K       '         Trip)  EE^n 

where  the  trace  is  over  both  orbital  indices  //,  v  and  spin  index  7,  and  S^  is  the 
overlap  matrix  element  in  this  basis.  The  denominator 


is  the  total  electronic  population,  and 


(3.4) 


"7  =  E  P7"S' 


(3.5) 


HV 


is  the  average  electronic  population  for  spin  7. 

The  electric  multipoles  contain  the  information  about  the  spacial  distribution 
of  the  electron  cloud.  For  k=Q,  p(°)  =  1.  For  k  >  0, 


is  the  component  of  the  dipole  moment  P  (since  the  electron  mass  is  1  au)  and 

(3.7) 


(2)_Tr[^(r2^-3r,rJ)]        1 

Tt\p\ =  n  L>  1"  P^T  8il  ~  3r,r^ 


PK  '  = 

*3 


VH 


7      *»" 


is  the  component  of  the  quadrupole  tensor  j  of  the  electron  cloud. 

The  electric  multipoles  can  be  defined  either  in  the  collision  frame  or  in  the 
natural  frame.  It  is  also  useful  to  define  them  in  a  body-fixed  "atomic  frame" 
centered  at  the  scattered  atom,  in  which  the  axes  are  along  the  major  axes  of  the 
electron  cloud  and  the  quadrupole  matrix  appears  diagonal. 
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3.3  Alignment  and  Orientation  Parameters 


Another  way  to  analyze  the  orbital  polarization  is  to  examine  the  anisotropy 
of  the  distribution  of  the  orbital  angular  momentum  of  the  electron  cloud.  For  the 
orbitals  of  certain  total  angular  momentum  J,  if  the  distributions  of  the  orbitals  for 
all  \JM)  states  are  the  same  the  cloud  is  isotropic;  if  the  orbital  distributions  on 
these  states  are  different,  the  cloud  is  said  to  be  oriented;  if  the  orbital  distribution 
for  \JM)  and  \J  -  M)  are  equal  but  differ  for  different  magnetic  number  M  the 
cloud  is  said  to  be  aligned.  To  precisely  describe  these  polarization  properties, 
we  define  the  orientation  and  alignment  parameters  through  rotational  moments 
[Fano  and  Macek,  73;  Blum,  81]. 

Introduce  a  set  of  spherical  irreducible  tensor  operators 


M  _  i 

Jq       -  1 

J±1  "  T^J± 

J(1)  -  J 

7(2)  _   1  /2 

J±2  —  2   ± 

41  =  T\j±(2Jz 

±d 

^  =  >- 

j2) 

(3.8) 


where  /  is  the  angular  momentum  operator  and  Jz  is  its  z-component,  /+  and  J_ 


are  the  raising  and  lowing  operators.  The  irreducible  tensor  operator  j\  '  satisfies 
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the  following  relations  [Zare,  88] 


-Mr 


(*) 


=  qJi 


(h) 


(3.9) 


J±,J, 


(*) 


h(*) 


=  [fc(*  +  1)  -  9(9  +  1)]»  JJS 


J. 


2]' 


(3.10) 


(3.11) 


(i) 


We  define  the  orientation  vector  0\  '  and  the  alignment  tensor  A\ '  as 


,(2) 


O?  V)  - 


4V) 


i 


>/j(  j + 1) 


i2e(J,(1)) 


(3.12) 


V6 


J8e(J}7) 


(2)v 


(3.13) 


v/J(J  +  l)' 
where  <  >  indicates  the  ensemble  average  for  a  given  /.  Using  the  density  operator 

and  density  matrix,  the  average  of  the  irreducible  operator  can  be  written  as 


(J?>>  = 


Tr(p) 


~EE^ 


J") 


7      H» 


Ufl 


(3.14) 


>(l) 


The  vector  0\  '(J)  describes  a  preferred  direction  of  angular  momentum.  The 

(2) 

tensor  A\"'{J)  describes  the  preferred  spatial  distribution  of  angular  momentum 
corresponding  to  each  \M\.  Oft' (J)  is  also  called  the  orientation  parameter.  In  the 
natural  frame,  it  is  the  measure  of  the  average  component  of  angular  momentum 
perpendicular  to  the  collision  plane,  J±. 
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The  meaning  of  these  parameters  is  more  transparent  and  the  calculation  is 
easier  in  the  natural  frame.  But  some  collision  problems  are  more  conveniently 
dealt  with  in  the  collision  frame.  One  can  first  construct  the  electron  density  matrix 
in  the  collision  frame  and  then  calculate  orientation  and  alignment  parameters  in 
the  natural  frame  after  a  frame  rotation. 

The  density  matrix  contains  all  the  information  on  the  electron  states  and  can 
be  used  to  define  other  parameter  to  describe  the  orbital  polarization.  For  example, 
we  can  define  the  alignment  angle  a(J)  which  measures  the  angle  between  the 
major  axis  of  the  /  component  of  the  electron  cloud  and  the  Xn  axis.  Choosing 
the  basis  set  such  that  y.  =  JM,  we  define,  for  7=1  [Nielsen  and  Andersen,  87] 


1 


<*i\)  =  1>[*  +  «rg(Piitf)] 


(3.15) 


and  for  7=2 


1 


"(2)  =  «  [*  +  arg(P22,2o  +  ^20,22)] 


(3.16) 


where  M  =  -M.  We  can  also  define,  for  7>2,  the  octupolar  angle  shift  tj (J) 
For  instance,  for  7=2 


7(2)  =  \arg{P22t2i)  -  a(2) 


(3.17) 


The  parameters  L±,  a(J)  and  tj(J)  have  practical  meanings  and  can  be 
compared  with  experimental  date  to  test  the  theoretical  model  and  methods  used 
in  obtaining  the  density  matrix  [Panev  et.  al.,  87;  Andersen  et.  al.,  86]. 
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3.4  Multipoles  and  Alignment  and  Orientation  Parameters  in  a  Subsystem 


The  parameters  defined  above  describe  the  orbital  polarization  of  the  whole 
system.  Sometimes  we  are  more  interested  in  the  orbital  polarization  of  a 
subsystem.  For  example,  in  atom-surface  collision  experiments,  spectroscopies 
detect  scattered  atoms  at  the  distances  where  the  influence  of  the  surface  is 
negligible  and  only  the  polarization  of  the  atomic  orbitals  is  of  interest.  Even 
in  the  theoretical  description  of  the  collision  states  at  short  distance,  only  a  small 
part  of  the  surface  should  be  taken  into  account,  since  the  atom  usually  interacts 
mainly  with  a  small  region  of  the  surface.  Let's  divide  a  system  into  a  primary 
region  and  a  secondary  region  which  are  labeled  by  p  and  s  respectively.  In  the 
case  of  atom-surface  collisions,  the  primary  region  can  be  chosen  to  contain  the 
scattered  atom  and  a  small  part  of  the  surface  and  the  interaction  between  the 
atomic  states  and  those  in  the  secondary  region  is  assumed  to  be  small.  In  the 
natural  frame  the  ensemble  average  of  J^  for  the  whole  system  can  be  written  as 

=  DEE^(44)).,-hE£^(44)) 

7       HEP  vGp  w        fi£P  j/gs 


i//i 


V\i 


fi£s  v£s 


(3.18) 

The  last  term  comes  from  the  secondary  region  and  is  not  of  interest.  The  second 
and  third  terms  contains  the  mixing  between  the  primary  region  and  the  secondary 
region. 
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Because  we  are  interested  in  the  orbital  polarization  of  the  scattered  atom  we 
define  a  set  of  irreducible  operators 

40)  =  i 

L{1)  -  L 

r(2)_   1r2 
L±2  ~  oL± 

L±l  =  fl:L±(2Lz  ±  1) 


(3.19) 


(2) 


V6 


m  -  l2) 


: 


where  L  is  the  orbital  angular  momentum  of  the  scattered  atom,  L?  its  z- 
component,  and  L+  and  L_    the  raising  and  lowing  operators.    The  average  of 


4k)  is 


(4.)  =  !^  =  J_EEE^(4,) 


Tta(p) 


1    H&A  v€p 

where  trace  is  over  the  subspace  of  atomic  states,  and 


vn 


(3.20) 


>vn 


(3.21) 


nA  =  TvA(p)  =  E  E  E  Pl*S» 

t    n€A  v&p 

is  the  average  electron  population  on  the  atom.  The  contribution  to  the  electron 
cloud  associated  with  the  secondary  region  is  not  included  since  the  mixing 
between  the  second  region  orbitals  and  the  atomic  states  is  small.  In  this  way  we 
can  define  the  orientation  and  alignment  parameters  for  the  scattered  atom  as 


0?\L)A  = 


1 


y/L(L  +  l) 


Re(L[1]) 


(3.22) 


40 


!A  y/L(L+l) 


i?e(42)) 


(3.23) 


The  alignment  angle  and  octupolar  angle  shift  can  also  be  defined  in  the  similar 
way.  At  small  distance,  the  surface  electrons  contribute  to  the  atomic  orbital 
polarization.  At  large  distance,  mixing  between  the  atomic  orbitals  and  the  surface 
orbitals  becomes  small  and  the  summation  of  v  vanishes  unless  for  v  e  A,  the 
polarization  parameters  are  equal  to  those  defined  for  an  isolated  atom. 


CHAPTER  4 
LINEARIZATION  OF  TDHF  EQUATIONS 


The  TDHF  equation  developed  in  Chapter  2  can  be  solved  straightforwardly 
only  for  simple  systems  requiring  small  basis  sets.  If  the  system  is  complicated, 
the  TDHF  equations  can  instead  be  solved  by  some  linearization  procedures. 
However,  one  has  to  be  very  careful  when  developing  a  linearization  method  for 
cases  where  the  density  matrix  oscillates  rapidly  with  time,  such  as  in  collisions 
accompanied  by  a  charge  transfer  through  a  near-resonance  process,  since  the 
solutions  of  the  linearized  TDHF  may  converge  slowly  or  not  converge  at  all. 

In  this  chapter,  we  describe  a  local  time  linearization  procedure  for  solving 
the  TDHF  equations  for  the  near-resonance  charge  transfer  process.    It  tackles 
the  rapid  variation  of  the  density  matrix  by  defining  a  time-dependent  reference 
density  matrix  P°7(,).    In  section  4.1,  we  develop  the  local  time  linearization 
procedure  for  a  general  case  and  obtain  a  linearized  equation  for  P°1(t),  which  is 
a  solution  for  P«  when  the  ion  position  is  fixed  in  space,  and  a  linearized  equation 
for  matrix  Q7(,),  which  is  the  first  order  change  of  the  density  matrix  due  to  the 
time-dependent  perturbation  caused  by  the  nuclear  motion.  The  formal  solutions 
for  P°7(r)  and  Q7(|)  are  obtained  with  use  of  an  exponential  transformation.  In 
section  4.2,  we  apply  this  procedure  to  a  system  in  which  the  electron-electron 
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Coulomb  interaction  is  ignored.    For  this  special  case,  the  coefficient  matrices 
in  the  equations  of  P°7(r)  and  QT(f)  are  time  independent  and  the  analytical 
solutions  can  be  obtained. 

4.1  The  Linearization  Procedure 

In  a  collision  which  involves  charge  transfer  through  a  near-resonance  process, 
electrons  jump  back  and  forth  among  several  close  energy  levels  and  this  leads  to 
a  rapid  oscillation  of  the  density  matrix.  In  the  contrast,  the  change  of  the  density 
matrix  under  the  time-dependent  perturbation  caused  by  the  motion  of  the  nuclei 
is  relatively  slow.  We  separate  the  two  time  scales  in  what  follows. 

For  a  small  time  interval  t0  <  t  <  t\,  we  define  a  matrix  P°T(r)  to  satisfy 


the  equation 


•P°7(0  =  So  lfJP*1(«J  -  P^OFfco 1 


(4.1) 


and  the  initial  condition  V^{t0)  =  P^(t0),  with  F^  =  F^(t0)  and  S0  =  S(*0). 
Hence  P°7(r)  can  be  interpreted  as  a  solution  of  the  TDHF  equation  when  the  ion 
nuclei  are  fixed  at  RA(t0).  The  evolution  P^(r)  oscillates  with  time  reflecting 
the  charge  transfer  among  the  energy  levels  with  frequencies  proportional  to  the 
inverses  of  the  energy  differences  of  energy  levels.  On  the  other  hand,  if  the  time 
interval  is  small,  the  change  of  the  density  matrix  due  to  the  nuclear  motion  is 
small.  To  get  a  linearized  equation  for  the  this  change,  we  define 


Qt(t)  =  P^)-P°7(*) 


(4.2) 


43 


z-1 


and 


From  (2.32)  we  have 


AS-i(0  =  S"1(<)-So- 


AF7(*)  =  F7(<)  -  F% 


:-l 


(4.3) 


(4.4) 


AF(<)  =  AH  +  AeGT(0  +  A„G7(*) 


(4.5) 


The  matrices  AeG  and  AnG  are  related  to  the  electron-electron  interaction  and 
are  given  by 


A«GZ„(<)  =  £^a(0(HM)<o  +  3TAW(H«A)<0 

/cA 


(4.6) 


and 

AnG^(i)  =  X^OACHM),  +  ^(OA^M),  (4.7) 

/cA 

where  _(^i/|«A)t   and  (/zi/||/cA)t  are  the  Coulomb  and  Coulomb  symmetrized 
2-electron  integrals  respectively  when  the  scattering  atomic  core  is  at  £a(q)> 
Inserting  Eqs.  (4.1)— (4.5)  into  Eq.  (2.31),  we  obtain  to  first  order  in  Q7(/), 
»Q*(t)  =  So  ^JQ^O  -  Q^OFj^o1 

+So1AH(*)P°7(*)  -  P°7(<)AH(OtS0-1 
+P°T(^)F2;tAS-1(0  -  AS-^OFZP0^)  (4.8) 

+S-1  AeG7(OP07(<)  -  P^r(*)A.jST(*)tSJT* 
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To  introduce  a  shorthand  notation,  we  use  a  double  index  notation  to  replace 
the  single  index  notation  in  the  following.  Let  upper  case  Roman  letters  refer 
to  double  indices,  then,  K  =  (fiis)  and  L=  («A).  With  this  notation,  matrices 
P0l\  Q7,  AeG7  and  AnG7  become  column  matrices  V0y,  Q? ',  AeQ  and  AnQ. 
Defining  two  square  double  index  matrices  X(i)  and  y(t)  such  that 


XKlfo)  =  ^/iU,KX(t)  =  (flv\\K\)t 


VkiW  =  Y^att)  =  (H«A), 

we  can  rewrite  Eqs.   (4.6)  and  (4.7)  as 

*4&M  =  EMuriGK*)  +  yQKLQl(t)} 


*»GJc(t)  =  J2l^KLV°LJ(t)  +  AyKLV°/(t)] 


or  in  the  matrix  form, 

Aeg7=x0Qi  +  y0Qi' 


Angy  =  ax  v01  +  Ay  v01' 

where  X0  =  X(t0),   y0  =  y(t0)  and 

AX  =  X(t)  -  X0 


(4.9) 


(4.10) 


(4.11) 


(4.12) 


(4.13) 


(4.14) 


(4.15) 
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Ay  =  y(t)  -  y0 


(4.16) 


Taking  advantage  of  using  the  double  index  notation,  the  linearized  TDHF 
equations  (4.1)  and  (4.8)  can  be  written  in  a  neat  form 

iV°-r  =A17V°'r  +  ^rr'poV  _  -pOfA^  _  vo-,' Ai't\  (4  17) 


and 


iQ?  =A1'fQ"'+A'n'Q7'  -  QMT7t  -  Q^A7^*  +  Vy  (4.18) 


where  A"1"1  and  Ayy  are  square  matrices  in  double  index  notation  or  four  index 


matrices  in  single  index  notation, 

*2l  =  *ll,Kx  -  £  5ok  (*£,  Ma*  +  *&**£)  (4-19) 


Cn 


AkL  -  A,xu,kX  -  2^  S0li(iY%r,^\Pr12 

and  the  column  matrix  V1  is  a  driving  term  with  elements 

V\  =  [S0-1(AH  +  AnG)p°T  -P°T(AHt  +  AnG^So"1 


(4.20) 


(4.21) 


+as-1f2:p07  -  p^fJas-1]^ 

If  we  collect  VQ*<  and  V0^  into  a  column  matrix  V°,  Q7  and  Qy'  into  a  column 
matrix  Q,  V"1  and  X>7'  into  a  column  matrix  X>  and  let 


Ayi     A™ y 


A=  i 


(4.22) 
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then  Eqs.    (4.17)  and  (4.18)  become 

iV°  =  AV°  -  V°A* 


and 


(4.23) 


iQ  =  AQ-  QA*  +  V 


(4.24) 


These  are  the  linearized  TDHF  equations  for  the  density  matrix  in  t0  <  t  <  h.  It 
should  be  noted  that  the  matrix  A  is  generally  time-dependent. 

To  solve  the  above  equations,  we  let 

t 


U(t,t0)=T[exp<  - 
U(t,t0)*=T[exp{i 


I  A{t')dt' 

to 

t 
J A{tfdt' 


] 


(4.25) 


] 


(4.26) 


where  T  [exp{f  •  •  •}]  denotes  a  time-ordered  exponential  expansion.  We  consider 
transformations 

V\t)=U{t,t0)VQv{t)U(t,t0)t 


(4.27) 


and 


Q(t)  =U{t,t0)Qv{t)U{t,t0)t 


(4.28) 


Replacing  Eq.   (4.27)  into  the  left  hand  side  of  Eq.   (4.23)  and  Eq.   (4.28)  into 
the  left  hand  side  of  Eq.    (4.24)  we  have 

iV°(t)  =  AV°  -  V°A^  +  iU(t,  t0)Vl{t)U(t,  i0)f  (4.29) 
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and 


iQ(t)  =  AQ-  QA*  +  iU(t,  t0)Qv(t)U(t, Q1 


(4.30) 


Comparing  them  with  the  right  hand  side  of  Eqs.  (4.23)  and  (4.24),  we  obtain 


iVUt)  =  0 


(4.31) 


and 


iQv(t)  =  iU{t,  U>)-xV(t)  \U(t,  <o)f 
Their  solutions  are 


-1 


(4.32) 


(4.33) 


and 


Qv(t)  =  Qv(t0)  -  i  /#(*'k]r1k(t^iVk)tr1<a'  (4-34) 


to 


Using  Eqs.  (4.27)  and  (4.29)  and  noticing  that  Q(t0)  =  0,  we  obtain  the  formal 
solution's  for  V°(t)  and  Q(t) 


V°(t)  =  U(t,t0)-1V°(t0)[u(t,t0)r 


and 


(4.35) 


Q(t)  =  -M(t,t0)  J 'u{t\t0yxv{t')[u{t',toyyldt'u{t,tJ 

to 

t 
=  -i /w(m')z>(*'Mm')W 


(4.36) 
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4.2  The  Case  without  Electron-Electron  Interaction 


In  this  section,  we  will  show  the  use  of  our  linearization  procedure  by 
applying  it  to  a  special  case  where  the  electron-electron  is  ignored.  In  this  case 
G  =  0  and  F  =  H.  Since  the  different  spins  are  then  uncorrelated,  we  return  to 
the  notation  involving  k,  X,  fi  and  v,  and  drop  the  spin  index  7  to  write  Eqs. 
(4.23)  and  (4.24)  as 


tP°(*)  =  WP°(*)-P°(*)Wt 


(4.37) 


and 


where 


iQ(t)  =  WQ(<)  -  Q(i)Wt  +  D'(0 


W  =  So  "Ho 


H0  =  H(*0) 


D'(i)  =  Sq1  AH(<)P°(0  -  P°(0AHt(t)So"1 
+P°(f)HoAS-1(0  -  AS-1(OH0P°(0 


(4.38) 


(4.39) 


(4.40) 


(4.41) 


Unlike  in  the  general  situation,  the  coefficient  matrices  W  and  W*  are  time- 
independent. 

To  solve  this  equation  let  us  consider  transformations 


P°(<)  =  exp[-iW(t  -  t0)}P°B(t)exp\iW\t  -  *0) 


(4.42) 
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and 


Q(«)  =  exp[-iW(t  -  t0)]Qj3(t)exp^iW\t  -  *„)]  (4.43) 

Inserting  them  into  Eqs.  (4.37)  and  (4.38)  respectively,  we  get  equations  for  Pd° 
and  QD 


«'P8>(<)  =  0 


(4.44) 


and 


*'Qd(0  =  exp[iW(t  -  *0)]D(*)exp  -iW*(t  -  *„)]  (4.45) 


The  latter  has  a  solution  of 

t 


Qd(0  =  -  *  f  exp[iW(t'  -  t0)]D'(t,yjexp  -iW^(t'  -  t0)]dt'        (4.46) 


since  Qd(*o)  =  0.  From  the  above  equations  we  can  write  the  formal  solutions 
for  P°  and  Q 


P°(0  =  exp[-iW(t  -  to)]P0(t0)exp\iW*(t  -  t0) 


(4.47) 


and 


t 
Q(0  =  ~«  /  ™P  -tW(t  -  t'^B'(t'yxpUwUt  -  t'}]dt'  (4.48) 


We  assume  that  the  matrix  W  can  be  diagonalized  by  a  linear  transformation, 
that  is,  there  exists  a  matrix  L  such  that, 


W  =  LwL 


-l 


(4.49) 
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where  w  is  a  diagonal  matrix.   Then 


hi 


E^uWMoHL-1) 


4 


where 


(4.50) 


t 
rk\Pi{t  t0)  =  Jexp[-i(wk  -  wf)(t  -  t')]DXp(t')dt'  (4.51) 


The  formalism  developed  in  this  section  can  be  directly  applied  to  systems 
with  a  few  electrons.  We  have  used  it  to  calculated  charge  transfer  in  H+H+ 
collisions  [Runge  et  al.,  90].  However,  it  must  be  modified  when  applied  to 
extended  systems.  In  Chapter  6  we  describe  a  partitioning  procedure  which 
simplifies  the  treatment  of  the  scattering  of  atoms  by  surfaces. 


CHAPTER  5 
PARTITION  OF  EXTENDED  SYSTEMS 


For  an  ion-surface  system  which  involve  a  great  number  of  electrons,  the 
TDHF  equations  obtained  in  Ch.  2  usually  require  a  huge  basis  set.  Some 
technique  for  truncating  the  system  as  well  as  the  basis  set  must  be  developed 
to  solve  the  equations. 

Partitioning  techniques  have  been  used  to  treat  various  kinds  of  extended 
systems  [Lowdin,  70;  Ying  77;  Williams  et  al,  82;  Kirtman  and  de  Melo  81; 
de  Melo  et  al.,  87;  McDowell,  82,  85].  The  main  idea  is  to  divide  the  system 
into  two  regions  and  to  employ  some  relatively  accurate  treatments  to  deal  with 
electronic  interactions  in  a  primary  region,  which  is  of  main  interest,  while  using 
certain  approximations  in  a  secondary  region.  The  difference  among  a  variety  of 
methods  is  in  the  treatments  of  the  secondary  region  and  of  the  interaction  between 
the  two  regions.  For  some  systems  the  interaction  between  the  two  regions  is 
small  and  the  secondary  region  can  be  treated  as  a  small  perturbation.  For  other 
systems  the  secondary  regions  are  often  approximated  by  some  simpler  systems, 
for  example,  in  studying  Kondo  effect  [Kondo,  69;  Heeger,  69;  Anderson,  61],  an 
approximation  is  to  treat  the  conduction  electrons  surrounding  isolated  magnetic 

ions  as  a  free  electron  sea. 

- 
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For  ion-surface  systems,  it  has  been  found  that,  during  collisions,  electron 
rearrangement  occurs  mainly  within  local  regions  of  surfaces  [Grimley  et  al.,  83; 
McDowell,  85].  Based  on  this  observation  we  design  a  partitioning  technique  to 
truncate  the  ion-surface  system  into  two  parts  and  treat  electronic  interactions  and 
electronic  state  evolution  in  these  parts  differently.  This  partitioning  technique  is 
suitable  to  deal  with  time-dependent  processes  in  extended  systems  with  a  strong 
local  coupling.  It  requires  a  set  of  localized  basis  functions,  but  does  not  depends 
on  the  concrete  form  of  the  functions.  All  the  formal  derivations  in  this  chapter 
are  done  with  localized  basis  functions  which  can  be  of  any  form. 

In  Section  5.1  we  describe  our  partitioning  method  and  derive  the  effective 
linearized  TDHF  equation  for  the  density  matrix  in  the  primary  region,  which 
includes  a  driving  term  containing  the  coupling  with  the  secondary  region.  In 
Section  5.2,  we  introduce  an  approximation  to  the  secondary  region  and  the 
couplings  between  the  two  regions,  which  enables  us  to  solve  the  effective  TDHF 
equations  in  the  primary  region.  We  apply  this  method  to  a  simple  case,  where 
the  electron-electron  interaction  is  ignored,  and  obtain  analytical  solutions  for  the 
density  matrix  in  the  primary  region. 

5.1  Partition  Procedure 

We  divide  an  ion-surface  system  into  two  regions:  (i)  the  primary  region  that 
consists  of  the  scattering  ion  and  a  small  impact  area  of  the  surface,  where  the 
perturbation  by  the  ion  is  strongly  felt  during  the  collision;  (ii)  the  secondary 
region  which  consists  of  the  remainder  of  the  surface,  where  the  influence  of  the 
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ion  is  assumed  to  be  felt  indirectly  through  the  coupling  between  the  electrons  in 
the  two  pans  of  the  surface.  Considering  a  set  of  localized  basis  functions  for 
the  solid  and  the  ion,  and  letting  p  be  the  index  of  the  primary  region  and  s  the 
index  of  the  secondary  region,  the  matrix  P°  is  written  as 

/pO         pO    \ 

p0=Up  p!')  4 

and  other  matrices  are  partitioned  in  a  similar  way.  In  the  case  where  the  electron- 
electron  interaction  is  ignored,  Eqs.  (4.37)  and  (4.38)  are  then  split  into  eight  equa- 
tions for  Pjp(0,  P°ps(0,  P°ap(t),  P°ss(t),  Qpp(t),  Qps(*),  Qsp(0  and  Q88(<) 
iP°pp(t)  =  WppP°pp(t)  -  ppp(*)wp 

+wp,p!p(<)  -  pps(*)wstp 

iP°ps(t)  =  WppP%(t)  -  Ppp(*)Wj 

+Wp.Pi(<)  -  Pps(0wt 

ip°sp(t)  =  wspp°pp(t)  -  p;p(owj 
+w»p;p(o  -  pi(owjp 


pp 


(5.2) 


ps 


ss 


pp 


(5.3) 


(5.4) 


iP°ss(t)  =  WspP°ps(t)  -  P°ap(t)Wl 
+WS3P°as(t)  -  P°a3(t)Wl 
iQpp(t)=WppQpp(t)  -  Qpp^Wjp 

+WpsQsp(0  -  Qps(i)Wstp  +  Dpp(0 
iQPs(0=WppQp9(<)  -  QPP(/)WPS 

+WpsQss(i)  -  Qps(i)W|s  +  Dps(0 


--%psv/—  "  pp-^psv-/       ^ppv;  "  ps 


(5.5) 


(5.6) 


(5.7) 
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PP 


(5.8) 


(5-9) 


'pp 


pp 


'sp 


«'Q.p(0=W8pQpp(0  -  Qsp(*)wP 

+WssQsp(0  -  Qss(OWstp  +  D'sp(t) 

iQss(t)=WspQps(t)  -  Qsp(OWps 

+W8SQss(i)  -  Qss(t)Wl  +  B'ss(t) 
where  Dpp(i),  Dps(f),  B'sp(t)  and  B'ss(t)   are  the  submatrices  of  D  (r)  defined  in 
Eq.  (4.21).  For  example,  the  matrix  Dpp(£)  is  given  by 

DppW  =  (So1)ppAHpp(0P°pp(0  -PPp(0AHPp(0(So"1) 
+  (So1)ppAHps(0P°p(0-PPS(0AHtp(^)(So-1) 
+  (So1)/Hsp(i)P;p(i)  -  PPp(0AHps(i)(So-1) 
+  (So  1)psAHss(0Ps°p(0  -  P°ps(*)AHsts(*)(So  *) 
+PPP«HSPp(AS-1)pp(0  _  (AS-1)pp(0HOppPPP(0 
+Ppp(^H;ps(AS-1)sp(0  -  (AS-1)ps(i)HOspPPP(0 
+P°psWHLP(AS-1)pp(0  -  (AS-1)pp(0HOpsP°p(0 
+Pps(^HL(AS-1)sp(0  -  (AS-iyOHossPV) 

The  primary  region  where  charge  transfer,  energy  transfer  and  other  collision 
phenomena  take  place  is  our  main  interest.  But,  solving  the  effective  equations  in 
the  primary  region,  Eqs.  (5.2)  and  (5.6),  is  not  any  easier  than  solving  the  original 
linearized  TDHF  equations,  since  they  contain  the  density  matrix  elements  in 
the  secondary  region  and  have  to  be  solved  simultaneously  with  the  equations 
for  density  matrix  elements  involved  the  secondary  region.    Obviously,  some 


(5.10) 
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approximation  has  to  be  made  in  the  secondary  region  to  make  this  partition 

procedure  practical. 

5.2  The  Approximation  in  the  Secondary  Region 

When  considering  the  secondary  region,  we  notice  that  during  the  collision, 
the  charge  exchange  is  mainly  restricted  to  the  primary  region.  To  understand  it, 
let  us  consider  the  interaction  and  the  evolution  of  the  electronic  states  in  the  two 
parts  of  the  surfaces,  i.e.,  a  small  impact  area  on  the  surface  which  is  close  to  the 
scattering  ion  during  the  collision  and  the  remainder  of  the  surface.  During  the 
collision,  the  impact  area  experiences  a  strong  and  time-dependent  perturbation 
from  the  ion  and  the  evolution  of  the  electronic  states  in  the  area  is  significantly 
altered  from  that  of  the  unperturbed  surface.  As  the  result,  charge  in  the  impact 
area  may  move  to  or  from  the  ion.  In  the  remainder  of  the  surface,  the  direct 
interaction  with  the  ion  is  much  weaker,  the  electrons  in  this  area  feel  the  effect 
of  the  ion  indirectly  through  their  coupling  with  electrons  in  the  impact  area. 
Because  the  time  scale  of  charge  rearrangement  in  the  solid,  TrearT,  is  much  longer 
than  the  duration  of  the  collision  at  the  surface,  Toi,  at  the  collision  energies  of 
present  interest,  the  electrons  in  the  remainder  of  the  surface  do  not  have  time  to 
adjust  to  the  charge  rearrangement  in  the  impact  area,  and  the  electronic  states 
follow  approximately  the  evolution  pattern  of  an  unperturbed  system.  Thus,  it  is 
reasonable  to  assume  that  the  secondary  region  is  basically  unperturbed  by  the 
ion  during  the  collision.  In  this  approximation,  the  submatrices  of  Q  associated 
with  the  secondary  region,  Qps(t),  Qsp(t)  and  Qss(t),  are  found  to  be  negligible 
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compared  to  Qpp,  and  the  submatrices  P°ps(t),  P°sp(t)  and  P°ss(t)  can  be  replaced 
by  those  of  the  uncoupled  system,  Pps(t),  Psp{t)  and  Pss(t).  Other  submatrices 
are  treated  in  the  similar  way. 

For  an  unperturbed  surface,  the  coefficient  of  the  molecular  orbital  has  the 
form  of 


Ci7»(0  =  ci»(Un)exp[-iei(t  -  tin)] 
where  £,  is  the  energy  of  the  z'th  orbital.  From  Eq.  (2.11) 

Pltuit)  =    2*   Citl(t)civ{t)  =    ]P   %(4»)3»(*fa) 
i=occ  i=occ 

=  Piii/(Un) 


(5.11) 


(5.12) 


As  the  result  of  these  approximations,  we  need  to  solve  only  the  effective 
equations  in  the  primary  region 


«pJp(0  =  Wppp;p(o  -  vUtywl 


(5.13) 


and 


*"Qpp(<)=WppQpp(0  -  QPP(0Wtp  +  &   (t) 


(5.14) 


The  equation  for  Qpp(t)    is  coupled  with  the  secondary  region  through  terms 
within  B'pp(t)    which  is 

DppW  -  (So1)ppAHpp(0Ppp(0  -  P;p(t)AHtp(<)(S-i) 

+  PPP(0<P(AS-1)pp(0  -  (AS-^^COHoppPjpft) 
+Pp.('.-«)H{8P(AS-1)     (t)  -  (AS"1) (t)H0psP%(ttn) 


'pp 


(5.15) 
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Given  the  matrix  elements  S^1  and  H0,  the  driving  term  T)'pp(t)  can  be 
calculated  and  the  density  matrix  in  the  primary  region  can  be  obtained  either 
analytically  or  numerically  by  using  Eqs.  (4.50)  and  (4.51). 

It  should  be  pointed  out  that  the  validity  of  our  approximation  in  the  secondary 
region  depends  on  two  factors,  the  ratio  of  rCoi/rreaiT  and  the  size  of  the  primary 
region.  For  metal  surfaces  and  colliding  atoms  with  their  kinetic  energy  above 
leV,  rcoi/rrearr  is  small  and  our  partitioning  procedure  should  apply.  As  for  the 
size  of  the  primary  region,  the  larger  the  primary  region,  the  more  accurate  the 
result.  But  increasing  its  size  will  significantly  enlarge  the  basis  set  in  the  primary 
region  and  increase  the  computing  time.  In  Ch.  9  we  will  discuss  the  effect  of 
the  primary  region  size  in  more  detail. 


CHAPTER  6 
ELECTRONIC  BASIS  FUNCTIONS  AND  MATRIX  ELEMENTS 


The  use  of  the  localized  basis  functions  in  the  present  study  appears  nec- 
essary for  our  partition  procedure  introduced  in  Chapter  4,  however,  the  real 
reason  underlying  it  is  the  local  nature  of  the  phenomena  associated  with  sur- 
faces. For  example,  localized  surface  states,  surface  adsorptions,  chemisorption 
and  atom-surface  collisions  have  their  effects  basically  confined  to  small  regions. 
Experimental  results  have  also  provided  evidence  that  charge  densities  and  local 
densities  of  states,  while  remarkably  different  from  those  of  bulks  on  the  top  lay- 
ers of  surfaces,  quickly  recover  to  the  bulk  values  on  deeper  layers.  For  example, 
the  electron  energy  distribution  results  obtained  from  the  ion  neutralization  spec- 
troscopy, which  probes  only  the  top  layer  of  atoms,  show  a  qualitative  difference 

from  those  of  the  bulk  [Hagstrum  and  Becker,  73],  but  the  results  from  the  ul- 

1 

traviolet  photoemission  spectroscopy,  which  probes  about  four  atom  layers,  are 
essentially  dominated  by  the  bulk  density  of  states  [Eastman  and  Grobman,  72]. 

Localized  basis  functions  have  been  used  in  calculating  electronic  structure 
of  bulk  materials;  the  examples  are  Wannier  functions  and  atomic-like  basis 
functions  used  in  the  tight-binding  method.  The  local  nature  of  the  surface 
phenomena  suggests  that  some  of  the  electronic  behaviors  of  surfaces  could 

58 


59 


be  more  advantageously  and  conveniently  described  in  terms  of  localized  basis 
functions  which  return  to  the  form  of  the  bulk  functions  on  deeper  layers. 

Recently,  some  works  have  reported  using  atomic  orbitals  to  calculate  surface 
electronic  properties.  Smith  and  his  colleagues  use  a  set  of  atomic-like  basis 
functions  to  calculate  surface  band  structure  for  transition  metals  [Smith  and 
Gay,  75;  Smith  et  al.,  80;  Alinghaus,  et  al.,  80].  Their  basis  functions  are 
constructed  by  fitting  them  to  the  solutions  found  by  solving  the  Kohn-Sham 
equations  [Kohn  and  Sham,  65].  Kohn  and  Onffory  introduce  the  concept  of 
the  generalized  Wannier  functions  which  can  be  used  in  the  theoretical  study 
and  calculation  of  electronic  structures  for  extended  systems  lacking  translational 
symmetry  [Kohn  and  Onffroy,  73].  These  functions  are  localized  on  the  lattice 
sites  and,  as  they  result  form  a  unitary  transformation  of  the  eigenfunctions  of 
the  system,  are  orthonormal  and  complete.  Gay  and  Smith  [Gay  and  Smith,  74] 
proposed  a  variation  method  for  determining  the  generalized  Wannier  functions 
without  having  to  first  calculate  the  eigenfunctions  of  systems. 

This  chapter  is  devoted  to  constructing  localized  functions  as  our  basis 
function  set  for  atom-surface  systems,  which  consists  of  a  set  of  localized  functions 
for  the  surface  and  a  set  of  atomic  functions  for  the  atom,  and  to  calculating  the 
relevant  matrix  elements  in  this  basis.  Although  the  localized  basis  functions  used 
by  Smith  and  his  colleagues  are  successful  in  giving  good  results  for  transition 
metals,  their  calculation  and  application  are  complicated.  As  a  test  of  our  time- 
dependent  molecular  orbital  method  and  the  partition  procedure,  we  feel  that  a  set 
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of  simple  localized  functions  would  be  adequate  and  hence  choose  the  generalized 
Wannier  functions  as  a  part  of  our  basis  set. 

In  Section  6.1,  the  definition  of  the  generalized  Wannier  functions  is  intro- 
duced and  a  variation  procedure  for  determination  of  these  functions  is  described. 
In  Section  6.2,  we  apply  the  variation  procedure  to  a  jellium  surface  to  obtain  the 
generalized  Wannier  functions.  For  convenience  in  the  calculation  of  the  matrix 
elements,  the  generalized  Wannier  functions  are  written  as  linear  combinations  of 
Gaussians.  In  Section  6.3,  the  atomic  functions  for  the  atom  are  constructed  by 
using  a  pseudopotential  for  the  Sodium  atom  core.  In  Section  6.4,  we  calculate 
the  matrix  elements  of  the  overlap  and  Hamiltonian.  Some  of  the  details  of  the 
calculation  are  given  in  Appendices  A  and  B. 

6.1  Generalized  Wannier  Functions 
6.1.1  Definition  of  Generalized  Wannier  Functions  (GWFs) 

We  first  consider  an  infinite  periodic  system  with  a  Hamiltonian  H°.  Its 
eigenfunction  <f>°Ar)  is  a  Bloch  function  labeled  by  the  wave  vector  k  and  band 
index  i  and  satisfies  the  Schrodinger  equation 


WW  =  wwn 


(6.1) 


where  e0-.  is  the  eigenvalue  associated  with  /  and  k.   The  functions  </>°-(f)  are 
then  collected  in  a  row  matrix 


$ 


'-(**»*A»"  '  •) 


(6.2) 
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and  the  above  Schrodinger  equation  has  a  matrix  form 


rO^O 


iOTTiO 


ff>$»  =  $YEY 


(6.3) 


where  E?  is  a  diagonal  matrix  with  diagonal  elements  e0-. 

ik 

For  each  energy  band  of  a  periodic  lattice,  Wannier  functions  are  defined  by 
[Wannier,  37] 


1 


"*«  -  m  E  •"*V^ 


\/JV 


(6.4) 


tg5Z 


where  N  is  the  number  of  the  lattice  points,  Rg  is  the  lattice  vector  labeled 
by  a  vector  index  n,  the  summation  is  over  the  Brillouin  zone.  The  Wannier 
function  w^(f)  is  localized  about  the  lattice  point  Rg,  In  matrix  form,  the  above 
transformation  can  be  written  as 


0    TTOt 


w,°  =  $JU 


where 


,o 


and 


«*-(«fc.«fc.  •  •  •) 


uhm'JJ[§**fi'**) 


(6.5) 


(6.6) 


(6.7) 


Since  matrix  U°  is  unitary,  the  Wannier  functions  are  orthonormal  and  are  an 
alternative  set  of  basis  functions  in  electronic  structure  calculation  [See,  for 
example,  Slater  and  Koster,  54]. 
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For  a  semi-infinite  solid  with  its  surface  parallel  to  the  x-y  plane,  its  Hamil- 
tonian  H  and  eigenfunction  (f>^(r),  which  is  characterized  by  the  band  index  i 
and  the  wave  vector  k,  satisfy  the  Schrodinger  equation 


**a(*0  =  'i&aW 


(6.8) 


where  e^  is  the  eigenvalue.  Collecting  the  functions  <t>&{r)  in  a  row  matrix 


>.=  (V.fc^.v-  •  •) 


(6.9) 


the  above  Schrodinger  equation  can  be  written  in  matrix  form 


##j  =  $,E, 


(6.10) 


where  E,  is  a  diagonal  matrix  with  diagonal  elements  tA  In  the  following,  we 
drop  the  band  index  i  with  an  understanding  that  all  the  functions  and  matrices 
involved  are  referred  to  a  particular  band. 

As  a  results  of  breaking  translation  symmetry  in  the  direction  perpendicular 
to  the  surface,  the  eigenfunctions  <^(r)  are  no  longer  of  Bloch-type  and  the 
functions  defined  in  Eq.  (6.4)  are  no  longer  meaningful.  Kohn  and  Onffroy  point 
out  that  it  is  possible  to  construct  a  set  of  orthonormal  functions  for  systems  with 
defects,  called  generalized  Wannier  functions  [Kohn  and  Onffroy,  73].  Following 
their  idea,  we  define  the  generalized  Wannier  functions  of  an  electronic  band  as 
a  unitary  transformation  of  the  eigenfunctions  of  the  system,  that  is, 


w  =  $Ut 


(6.11) 
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where 


W  =  (Wft^Wfta,   ■    •    •) 


(6.12) 


is  the  row  matrix  of  generalized  Wannier  functions  and  U  is  a  unitary  matrix  to 
be  determined.  Obviously,  w  is  orthonormal,  i.e., 


(w|w)  =  I 


(6.13) 


It  has  been  proved  that  Generalized  Wannier  functions  are  localized  about 
the  lattice  sites  fig  and  are  exponentially  decaying  away  from  their  center  [Kohn 
and  Onffroy,  73].  The  well-behaved  localization  property  of  the  generalized 
Wannier  functions  combined  with  their  orthonormality  make  them  particular  useful 
to  investigate  the  electronic  structures  or  localized  phenomena  for  non-periodic 
systems. 

6.1.2  Generalized  Wannier  Functions  as  Linear  Combinations  of  Gaussians 

For  the  purpose  of  practical  applications  we  write  the  generalized  Wannier 
functions  as  linear  combinations  of  Gaussian  functions  centered  at  different  sites. 
For  a  certain  energy  band  only  Gaussian  functions  with  proper  symmetry  should 
be  used  to  construct  the  generalized  Wannier  function.  Thus, 


m 

where  Bn^  is  a  coefficient,  a=(nlm)  is  a  compound  index,  and 
<7am(/?am,r-)  =  baYtm(9,  ip)rl  e~  ^^  ^ 


(6.14) 


(6.15) 
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is  the  a  type  Gaussian  primitive  function  centered  at  R^,  where  F/m(0,  tp)  is  the 


spherical  harmonic  functions  and  ba  a  normalization  factor.  Letting 

o         \9arh\  i   9am2 1    '    '    '  ) 


(6.16) 


(here  we  have  dropped  index  a  in  g  since  a  is  the  same  for  all  the  Gaussians  for 
a  given  band),  Eq.  (6.14)  can  be  written  in  a  matrix  form 


w  =  Bg 


(6.17) 


where 


(B)nm  =  Bnm 


(6.18) 


Substituting  Eq.  (6.17)  into  the  normalization  relation  Eq.  (6.13)  yields 


where 


BfGB  =  I 


G  -  (g|g) 


(6.19) 


(6.20) 


is  the  overlap  matrix  of  Gaussians,  we  find  that  a  possible  choice  for  B  is 


B  =  G" 


(6.21) 


where  the  matrix  G~2  satisfies 


G-*GG-*  =  I 


(6.22) 
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6.1.3  Determination  of  Generalized  Wannier  Functions 

The  generalized  Wannier  functions  can  be  constructed  from  the  eigenfunctions 
<f>%(r)  according  to  their  definition  (6.11).  However,  one  of  the  important  features 
of  the  generalized  Wannier  functions  is  that  they  can  be  directly  determined  from 
the  system  Hamiltonian  without  having  to  know  the  original  eigenfunctions  [Kohn 
and  Onffroy,  73;  Gay  and  Smith,  74],  which  allows  freedom  to  construct  them 
as  desired.  In  this  study  we  use  a  variational  procedure  [Gay  and  Smith,  74]  to 
find  the  generalized  Wannier  functions  w. 

Since 

E»{#|&|#) 


(6.23) 


(6.24) 


the  total  energy  of  a  band 

EB  =  Tr(E)=Tr((&\H\$)} 
=  7/r(<w|#|w))=£fl[w] 

is  a  functional  of  the  w's,  and  we  have  used  in  the  second  line  the  orthonormality 
Eq.  (6.13).  The  variational  principle  proposed  by  Kohn  and  Onffory  [Kohn  and 
Onffroy,  73]  says  that  the  total  energy  of  a  band  attains  its  minimum  if  a  correct 
set  of  w's  is  used.  Thus  the  generalized  Wannier  functions  can  be  determined  by 
minimizing  £fl[w],  that  is, 


min{£5[w]}  =  £ 

{w} 


while  being  subject  to  the  constraint  Eq.   (6.13). 


(6.25) 
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Once  w  is  known  we  can  insert  Eq.    (6.11)  into  the  eigenequation  of  the 


system  to  have 


#wU  =  wUE  (6.26) 


or,  using  Eq.    (6.13), 


(w|tf  |w)U  =  UE  (6.27) 


This  equation  can  be  used  to  find  the  matrix  U  as  well  as  the  eigenfunction  <f>r(r) 


By  writing  the  generalized  Wannier  functions  as  linear  combinations  of 
Gaussians,  the  functional  EB[w]  becomes  a  function  of  fa%  n(0h  fc,  •  •  •)> 
and  the  problem  of  searching  for  a  set  of  w  becomes  one  of  finding  a  set  of 


optimized  0**9,  i.e., 


min    [n(ft,  ft,  •••)]-»  A,  fo,  ■  •  ■  (6.28) 

Pl,P2,— 


One  may  recall  that  the  above  procedure  for  writing  the  generalized  Wannier 
functions  as  linear  combinations  of  Gaussians  is  similar  to  those  used  in  molecular 
orbital  calculation.  However  the  difference  is  that  the  Gaussian  function  or  the 
exponent  fa  depend  on  the  site  position  R^.  For  a  surface,  due  to  the  periodicity 
in  x  and  y  directions,  all  the  generalized  Wannier  functions  on  the  same  layer  are 
the  same  and  the  exponent  fa  changes  only  with  layers.  From  now  on  we  will 
denote  fa  by  flmz ,  where  mz  is  the  z  component  of  the  vector  index  m,  to  indicate 
that  it  is  the  exponent  for  the  mzth  layer. 
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Before  proceeding  to  find  0mx's  we  notice  an  asymptotic  behavior  of  the 
generalized  Wannier  functions;  they  approach  the  Wannier  functions  of  a  periodic 
system  in  an  exponential  manner,  as  proved  by  Kohn  and  Onffory  [Kohn  and 
Onffroy,  73].  This  indicates  that  only  a  small  number  of  the  generalized  Wannier 
functions  need  to  be  determined  for  near-periodic  lattices  which  loose  translation 


symmetry  only  in  local  areas  or  in  certain  directions,  such  as  imperfect  crystals  and 
surfaces.  For  a  perfect  bulk,  the  Wannier  functions  w^(r)  's  on  all  the  sites  R&  are 
the  same  and  have  the  same  exponent  parameter  0°.  For  a  semi-infinite  lattice 
with  its  surface  at  m2=l,  the  generalized  Wannier  function  w^(f)  approaches 
w%(f)'s  and  0m,  approaches  0°  when  R&  is  deep  inside  the  surface.  There 
exists  a  cutoff  M  such  that  when  mz  >  M  it  is  found  approximately  that 

An.  =  P°  (6.29) 

The  variation  procedure  (6.28)  now  becomes 

min    [()(&,  #>,...)]=      min      [0(01,  ft,  ■  •  ;M]  -  A.A>  •  •  -t0M   (6.30) 
m»»i—  ri,P2,—pm 

Thus  only  M  variational  parameters,  0U  fa,  •  •  •,  0m,  must  be  determined. 
The  choice  of  the  number  M  depends  on  the  system  and  the  accuracy  requirement 
of  calculations. 

6.2  Generalized  Wannier  Functions  for  a  Jellium  Surface 

The  method  developed  in  the  previous  section  can  be  applied  to  construct 
generalized  Wannier  functions  for  any  system  for  which  the  effective  one  electron 
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Hamiltonian  is  known,  although  the  calculation  of  the  three  dimensional  matrix 
G~2  is  a  little  time-consuming  [Gay  and  smith,  74].  In  this  section  we  apply 
the  method  to  a  simple  model,  a  jellium  slab  with  a  step  potential.  A  jellium 
is  a  system  with  the  atom  core  lattice  replaced  by  an  uniform  positive  charge 
background  and,  as  a  result,  it  completely  ignores  the  effect  of  the  atomic  lattice. 
Because  of  its  simplicity  in  the  electronic  structure  calculation,  jellium  is  often 
used  in  condensed  matter  physics  as  a  testing  model  for  theory  or  methods.  We 
also  assume  that  there  is  only  one  band  and  use  Is  type  Gaussians  to  construct 
generalized  Wannier  functions.  We  fit  the  band  parameters  such  as  Fermi  level, 
and  the  energy  of  the  bottom  of  the  band  to  those  of  W(110). 

For  a  finite  jellium  slab  with  a  thickness  D  in  z-direction  and  a  step  potential 

(6.31) 


\     ° 

z>0 

V(r)  = 

-v. 

-D  <  z  <  0 

o 

z  <  -D 

the  Hamiltonian  is 


H(f)  =  -1-V2  +  V(r) 


(6.32) 


2 

=  Hx(x)  +  H*(y)  +  Hz{z) 
where  Hx(x)  and  Hy(y)  equal  to  the  kinetic  energies  in  x  and  y  directions  and 


«"W--j|^  +  VW  (6.33) 

with 

Vz(z)  =  V(r)  (6.34) 

Labeling  x,  y  and  z  components  by  superscripts  x,  y  and  z,  the  eigenfunction  Mr), 
which  is  assumed  to  satisfy  cyclic  boundary  conditions  in  x-  and  y-directions,  and 
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eigenvalue  Er  have  the  forms  of 


W-*<*MLto*M 


(6.35) 


and 


Ej-  =  Ekx  +  Eky  +  Eki  (6.36) 

respectively.  In  Appendix  C  we  describe  the  evaluation  of  the  Fermi  energy. 

We  now  consider  an  artificial  cubic  mesh  in  the  slab,  which  contains 
Ns=NxNyNz  sites  with  a  distance  d  between  sites,  where  Nx,  Ny  and  N2  are  the 
numbers  of  sites  in  x-,  y-  and  z-directions,  respectively.  The  generalized  Wannier 
functions  for  the  slab  are  associated  with  sites  and  are  written  as 


«>n(0  =  <,(x)w'(y)w'nt(z) 


(6.37) 


Recalling  Eq.  (6. 1 1),  the  matrix  U  can  then  be  written  as  a  product  of  x,  y  and 
z  component  matrices 

U  =  UxUyUz  (6.38) 


with  Ux,  Uy  and  Uz  lead  to  the  transformations 


*rr*  v        /  kxnx 

kx 


<w-?MJw 


(6.39) 


(6.40) 


<w  =  E  (v«)k,JiM 


(6.41) 
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respectively.   For  this  system,  therefore,  a  three  dimensional  calculation  is  sim- 
plified to  a  one  dimensional  one. 

For  this  system,  the  x-  and  y-components  of  the  eigenfunctions  are 

-  1  la 

(6.42) 


<!>kx  = 


y/Wt  - 1)3 

the  generalized  Wannier  functions  in  the  x-  and  y-directions  are  just  the  Wannier 

function  for  periodic  lattice 


Kx(x)  = 


I             sin 

ir(x— nxd) 
d 

(Nx-l)Vd    at-n 

T(x—nxd) 
[  (Nx-l)d  J 

(6.43) 


sin 
n"KyJ      (Nx-l)Vd    sin 


Tr(y-nyd) 
a 


r(y-nyd) 
(Nx-l)d 

and  the  component  transformation  matrices  Ux  and  l^  are 


(6.44) 


Unxkx  ~ 


U»         = 


1 


y/N^T 


exp(ikxnxd) 


(6.45) 


1 


*»*»  =  JN  _lex^-ikvnv^ 


(6.46) 


In  the  z  direction,  the  generalized  Wannier  function  for  «zth  layer,  u%  (z),  is 
written  as  Unear  combination  of  Gaussians 


<(^)  =  EBU4(/5m^) 


(6.47) 


m. 


71 
where  the  normalized  one  dimensional  Gaussian  function  in  the  z  direction 
Sm„(A»*«»*)  *s  centered  at  Zmi  on  mzth  layer  with  an  exponent  f3mz,  and 
Bz  =  (Gz)  2  with  Gz  =  (gz|gz)  the  overlap  matrix  of  Gaussians  in  the  z  di- 
rection. As  discussed  in  the  previous  section,  only  the  exponent  parameters 
on  the  first  M  layers  are  assumed  to  be  different;  that  is,  /?i  =  /?#,,  fa  = 
Pnz-1,   ■    ■    -,0M-1  =  Pn,-M+2,  but  @m  =  0M-1  =    •    •    •    =  @N,-M+1- 

In  the  x  and  y  directions,  the  Wannier  functions  on  the  mzth  layer  are  written 
as  linear  combinations  of  one  dimensional  Gaussians  with  an  exponent  /?mj.  For 
example,  the  Wannier  function  in  the  x  direction  is 

<m'(*)  =  E  *»,m,(A».)£,(A»., «)  (6-48) 

where  the  superscript  mz  for  the  function  is  used  to  indicate  the  layer  the  function 
is  on.  The  coefficient  matrix  element  #£imi(/?mJ  is  obtained  by  the  least  square 
principle  [Shavitt,  63]  by  fitting  Eq.  (6.48)  with  Eq.  (6.43).  The  details  of  the 
calculation  are  given  in  Appendix  A. 

From  Eq.    (6.24)  the  total  energy  is 
Eb  =  X>«|HK) 

n 

=  E  E  E  B**B**>(9*(0m.M9*>(fim'.))  <6-49> 

n       m     m' 

with  the  details  of  the  calculation  given  in  Appendix  B.  Minimizing  Q,  respect  to 
the  /?'s  gives  a  set  of  optimized  /3's. 
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To  determine  the  matrix  LP,  we  substitute  Eq.  (6.41)  into  the  eigenequation 
of  4l(g) 


H'tf  ,(*)  =  EkJl(z)  (6.50) 


to  have 


£  Vk,m,E'w'm,(z)  =  Ek,  £  &*.*,<,(*)  (6.51) 


m, 


Multiplying  it  by  w^z(z)*  and  integrating  over  z  give  the  equation  for  transfor- 
mation matrix  Uz 

J>fc,m,«|H>U  =  Ek.i?k,nM  (6.52) 

6.2.1  Results 



We  have  calculated  the  generalized  Wannier  functions  for  a  finite  jellium 
slab  modeling  W(110),  with  the  depth  of  the  potential  well  V0=0-846  au.  The 
distance  between  sites,  d,  is  not  equal  to  the  lattice  constant  of  W,  and  is  chosen 
to  be  d=2.6l  au  so  that  the  average  number  of  electrons  in  the  volume  associated 
with  a  site  is  one.  The  optimum  exponents  of  the  Gaussians,  /3m2,  which  are 
determined  by  minimizing  the  total  energy  Eq.  (6.49),  are  listed  in  Table  6.1.  It 
has  been  found  that  the  optimum  exponent  decreases  for  the  first  couple  of  layers 
but  changes  very  little  after  the  third  layer.  This  shows  that  the  effect  of  the 
surface  on  the  generalized  Wannier  functions  is  basically  limited  to  the  several 
top  layers.   Gay  and  Smith  also  find  that  the  exponents  are  different  from  that 
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Table  6.1  The  optimum  Gaussian  exponents  for  the  first  three  layers. 


Layer 

First 

Second 

Third 

Pmz 

0.53 

0.52 

0.30 

of  the  bulk  only  on  the  top  two  or  three  layers  even  though  they  use  different 
potentials  [Gay  and  Smith,  75]. 

The  generalized  Wannier  functions  for  the  first  four  layers  are  plotted  in  Figs. 
6.1,  6.2,  6.3  and  6.4,  assuming  /?mz=0.30  for  mz  >  3.  These  figures  show  that 
the  generalized  Wannier  functions  are  well  localized  about  their  centers  with  long 
decaying  tails.  For  the  first  two  layers,  the  effect  of  the  surface  is  significant,  the 
generalized  Wannier  functions  are  asymmetric  and  their  tails  oscillate  little.  On 
moving  to  the  deeper  layers,  the  generalized  Wannier  functions  are  more  and  more 
characteristic  of  an  infinite  bulk.  For  the  fourth  layer  the  tail  of  the  generalized 
Wannier  function  towards  the  center  of  the  slab  is  very  similar  to  that  of  the  bulk 
Wannier  function.  Although  the  surface  is  a  severe  perturbation,  our  results  show 
that  the  generalized  Wannier  functions  rapidly  become  the  bulk  Wannier  functions 
when  moving  away  from  the  surface,  which  is  in  agreement  with  the  theoretical 
prediction  for  the  generalized  Wannier  functions  [Kohn  and  Onffroy,  73]. 

To  exam  the  validity  of  using  generalized  Wannier  functions  as  basis  functions, 
we  compare  the  eigenvalues  and  wavefunctions  calculated  from  these  functions 
with  the  exact  ones.  Table  6.2  lists  both  exact  and  calculated  eigenvalues,  Figs. 
6.5,  6.6  and  6.7  depict  both  exact  and  calculated  wavefunctions  for  Ekz =-0.8444, 
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-0.6906  and  -0.3990  for  states  jz=\,  10,  17  respectively.  The  agreements  are  good 
for  all  the  energies,  and  especially  for  lower  energies.  Since  the  point-by-point 
comparison  of  wavefuntions  is  probably  the  most  stringent  test  of  the  accuracy, 
the  good  agreement  in  our  calculation  will  assure  the  accuracy  in  the  calculation 
of  observables  such  as  charge  density. 


In  the  x  and  y  directions,  the  Wannier  functions  to*  (af)  and  Wny{y)  on  the 
mzth  layer  are  approximated  by  linear  combinations  of  Gaussians  with  the  mzth 
layer  exponent  /?mi .  The  accuracy  of  the  approximation  depends  on  the  number 
of  Gaussians  used.  It  is  found  that  for  the  larger  0mi  (close  to  the  surface  in 
our  case)  more  Gaussians  are  needed  to  achieve  the  same  accuracy.  Figs.  6.8, 
6.9  and  6. 10  show  the  exact  and  calculated  Wannier  functions  for  the  first  three 
layers  when  21  Gaussians  are  used. 

The  calculation  and  application  of  generalized  Wannier  functions  for  a  jellium 
slab  reveal  some  nice  features  of  these  functions.  Firstly,  they  are  well  localized 
about  their  centers  with  their  tails  characterizing  the  translation  symmetry  of  the 
systems.  The  functions  rapidly  recover  to  the  bulk  Wannier  functions  away 
from  defects  or  surfaces,  thus  only  a  few  generalized  Wannier  functions  need 
to  be  determined.  Secondly,  the  generalized  Wannier  functions  can  reproduce  the 
wavefunctions  and  eigenvalues  to  a  satisfactory  accuracy.  These  attractive  features 
make  generalized  Wannier  function  suitable  for  electronic  structure  calculations 
for  surfaces  or  other  systems  with  a  broken  translation  symmetry,  and  makes 
them  useful  in  studies  of  local  phenomena. 


75 


0.75 


V      0.50 

a 

i 

O    o.oo 


0.25  - 


-0.25 


1 P 


-12.0 


1 — p 


-9.0 


-i — p 


-i — p 


-6.0  -3.0 

z(D) 


0.0 


3.0 


Figure  6.1  Generalized  Wannier  function  wfu(z)  on  the  first  layer  of  a  jellium  slab 

with  a  square  well  potential  and  a  potential  depth  V0=0.846.  The  surfaces  of  the  slab 

are  located  at  z=0  and  z— 21d,  where  d=2.61  au  is  the  distance  between  layers. 
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Figure  6.2  Generalized  Wannier  function  w^t(z)  on  the 
second  layer  of  a  jellium  slab  with  a  square  well  potential. 
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Figure  6.3  Generalized  Wannier  function  w^  (z)  on  the 
third  layer  of  a  jellium  slab  with  a  square  well  potential. 
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Figure  6.4  Generalized  Wannier  function  w^(z)  on  the 
fourth  layer  of  a  jellium  slab  with  a  square  well  potential. 
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Table  6.2  Exact  eigenvalues  and  approximate  eigenvalues  calculated  by  using  generalized 
Wannier  functions  for  a  slab  with  a  square  well  potential  given  by  Eq.  (6.31). 


Exact  eigenvalue 

Approximate 
eigenvalue 

-0.8444 

-0.8444 

-0.8398 

-0.8397 

-0.8320 

-0.8320 

-0.8211 

-0.8209 

-0.8071 

-0.8068 

-0.7900 

-0.7895 

-0.7698 

-0.7690 

-0.7465 

-0.7453 

-0.7201 

-0.7185 

-0.6906 

-0.6883 

-0.6580 

-0.6549 

-0.6224 

-0.6181 

-0.5837 

-0.5778 

-0.5420 

-0.5340 

-0.4973 

-0.4864 

-0.4496 

-0.4349 

-0.3990 

-0.3792 

-0.3454 

-0.3185 

-0.2890 

-0.2517 

-0.2299 

-0.1769 

-0.1682 

-0.1299 

a 
o 

o 

g 


-25.0      -20.0      -15.0       -10.0       -5.0         0.0 

z(D) 


5.0 


Figure  6.5  Exact  and  approximate  wavefunctions  associated  with  the 
exact  eigenvalue  -0.8444  au  for  an  one-dimensional  square  well. 
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Figure  6.6  Exact  and  approximate  wavefunctions  associated  with  the 
exact  eigenvalue  -0.6906  au  for  an  one-dimensional  square  well. 
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Figure  6.7  Exact  and  approximate  wavefunctions  associated  with  the 
exact  eigenvalue  -0.3990  au  for  an  one-dimensional  square  well. 
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Figure  6.8  Exact  and  approximate  Wannier  functions  for  a  slab  in  the  directions  parallel 
to  the  surface.  The  latter  is  a  linear  combination  of  Gaussians  with  an  exponent  /?=0.53. 
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Figure  6.9  Exact  and  approximate  Wannier  functions  for  a  slab  in  the  directions  parallel 
to  the  surface.  The  latter  is  a  linear  combination  of  Gaussians  with  an  exponent  /3=0.52. 
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Figure  6.10  Exact  and  approximate  Wannier  functions  for  a  slab  in  the  directions  parallel 
to  the  surface.  The  latter  is  a  linear  combination  of  Gaussians  with  an  exponent  0=030. 
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6.3  Atomic  Basis  Functions 


Since  in  low  energy  atom-surface  collisions  only  the  valence  electrons  of 
atoms  are  actively  involved  in  charge  exchange,  it  is  a  good  approximation  to 
ignore  the  inner-shell  electrons  and  use  a  pseudopotential  for  the  atomic  core. 
In  this  section  we  construct  the  atomic  basis  functions  for  3s  and  3p  orbitals  of 
Sodium  and  a  pseudopotential  for  its  atomic  core. 

We  choose  Slater-type  orbitals  (STO)  <j)a  (r;  (,  RA)  as  basis  functions,  where 
a=3s,  3px,  3py,  and  3pz,  £  is  the  exponent  of  the  orbitals,  and  RA  is  the  position 
of  the  atomic  core.  For  the  convenience  of  calculation,  we  use  STO-KG  functions 
<j)a  If;  C,  Ra)  which  are  linear  combinations  of  K  Gaussians  and  satisfy 

4>*  (r;  (,  Ra)  =  (Ha  (f;  1,  RA)  (6.53) 

and  are  in  the  form  of 

^35 [r ;  1 ,  Ra)  =  J2  d^k  '  9u  (r ;  ak ,  RA^  (6.54) 


hPj  (f;  1,Ra)  =J2  d3ptk  ■  g2pj  (f;  ak ,  RA)  (6.55) 

where  j=x,  y  and  z,  and  gls  (f;  ak,  RAj  and  g2pj  (f;  ak,  rA  are  Is  and  2p  type 
Gaussian  primitive  functions  respectively, 

3      ~ 
gis(r;ak,RA"j  =  f-f^J    exp(-ak\f  -  RA\2) 


(6.56) 


92P(r;ak,RA)  =  (— -^*j    rjexp(-ak\f-  £a|2)  (6.57) 
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We  use  K=3  and  £=1.27.  The  data  for  ak,  d33jk  and  d3ptk  [Hehre  et  aL,  70]  are 
listed  in  Table  6.3. 


Table  6.3  Coefficients  and  exponents  for  Gaussians  in  STO-3G  functions. 


k 

"Jt 

d3s,k 

^3p,ifc 

1 

5.27266(-2) 

9.00398(-l) 

4.62001(-1) 

2 

1.34715(-1) 

2.25595(-l) 

5.95167(-1) 

3 

4.82854(-l) 

-2.19620(-1) 

1.05876(-2) 

6.4  Overlap  Matrix  Elements 

Our  localized  basis  set  consists  of  generalized  Wannier  functions  localized 
at  the  sites  of  the  slab  and  atomic  functions  centered  at  the  atomic  core  of  the 
colliding  ion  as  described  in  the  previous  sections,  that  is 


{<M  =  {{*»},    {*.}} 


(6.58) 


Since  both  {<^}  and  {<f>a}  are  orthonormal  sets,  the  diagonal  elements  of 
the  overlap  matrix  S  are  1  and  all  the  off-diagonal  elements  are  null  except 
San  y^Aj  =  (<f>a\4>n)  =  <S"na  f  R-A)  which  are  functions  of  the  atomic  core  position 
RA.   Using  Eqs.    (6.23),  (6.53)  and  (6.54) 

K 

San[RA)  =J2J2d^BnmGaih(ak^mz,RA)  (6.59) 

L 1        .a  ' 


*  =  1 
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where 


Garh  (o*,  fan,,  MA)  =  (ga  (a*,  RA)  Ja*(&0)  (6.60) 

is  the  overlap  matrix  of  Gaussians  and  is  also  a  functions  of  ^. 

6.5  Hamiltonian  and  Its  Matrix  Elements 

In  this  section  we  construct  the  one  electron  Hamiltonian  for  the  system  of  a 
Na  atom  and  a  jellium  slab  and  calculate  its  matrix  elements  in  the  basis  described 
in  the  previous  sections.  The  one-electron  Hamiltonian  of  the  system  is 

H  =  -^V2  +  VA(r,RA)  +  VM{r)  (6.61) 

where  VA(f, RAJ  is  the  atomic  potential  and  VM{r)  is  the  potential  of  the  slab 
which  we  choose  to  as  the  one  given  by  Eq.  (6.31).  However,  due  to  the  use  of 
the  localized  basis  functions  and  the  partition  method,  the  effective  Hamiltonian 
in  the  primary  region  should  be  reconstructed,  which  includes  the  construction 
of  the  pseudopotential  for  the  atomic  core  and  adding  a  correction  term  to  get 
correct  electronic  couplings. 
6.5. 1  Pseudopotential  for  the  atomic  core 

In  the  present  study,  since  we  are  not  intending  a  full  ab  initio  calculation, 
a  pseudopotential  for  the  atomic  core  is  used  for  VA  (f,  B.A  to  eliminate  the 
inner-shell  electrons  from  the  calculation  [Kahn  et  al.,  76].  The  pseudopotential 
usually  has  a  form  of 


/=o 


(6.62) 
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where 


m=—l 


(6.63) 


is  the  projection  operator  on  \lm)  and  the  functions  V1  can  be  chosen  in  one 
of  many  possible  forms  [Szasz,  85]  and  usually  contain  some  parameters  whose 
values  are  chosen  to  give  agreement  with  either  the  results  of  ab  initio  calculation 
or  with  that  of  experiments.  In  this  study  we  choose  the  form  suggested  by 
Schwartz  and  Switalski  [Schwartz  and  Switalski,  72], 


T///.  5   \  Z-Ne  exp(-K,\r  -  RA\A 

V'[r,RA)  =  -- ^-  +  At i ^ L  (6.64) 

V  J  r-RA  \r-RA\ 


where  Z  is  the  atomic  number,  Ne  the  number  of  electrons  of  the  core,  and  At  and 
ki  are  the  pseudopotential  parameters.  The  first  term  is  the  Coulomb  potential 
of  a  charge  of  Z-Ne  and  the  second  term  is  a  correction  due  to  the  electrons  in 
the  inner  shells. 

Using  the  pseudopotential  the  eigenequation  for  an  isolated  Na  atom  is 

"2 v2  +  va (r, RA)    <t>a  (r, RA)  =  ta<t>a (f, RA)  (6.65) 

We  determine  the  parameter  At  and  «/  by  fitting  the  eigenvalues  ea  with  the 
experimental  values  e3s=-0. 18884  au  and  e3p=-0.11156  au  [Callaway,  69].  An- 
other constrain  imposed  when  choosing  the  parameters  is  that  the  pseudopotentials 
should  become  Coulombic  at  the  distance  where  the  core  dies  off,  for  example, 
at  3r/,  where  3r/  is  the  radius  of  orbital  /.  However,  the  variation  of  parameter 
Ai  is  not  too  critical  [Schwartz  and  Switalski,  72]  and  can  be  chosen  in  a  cer- 
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tain  range.  The  parameters  for  our  pseudopotential  are  listed  in  Table  6.4.  The 
pseudopotentials  for  3s  and  3p  orbitals  are  plotted  in  Fig.  6.11. 

Table  6.4.  Pseudopotential  parameters  for  3s  and  3p  orbitals  of  Na  atom. 


Ai 

K[ 

3s 

1.0 

2.415 

3p 

3.0 

0.515 

6.5.2  A  Correction  Term  to  the  Hamiltonian 

The  partition  procedure  described  in  Ch.  5  can  be  used  to  describe  the 
evolution  of  the  electronic  states  and  the  collisional  charge  transfer  at  short 
distances  because  of  the  local  nature  of  the  problem.  But  it  complicates  the 
description  of  the  interaction  between  the  surface  and  the  atom  when  the  atom  is 
far  away  from  the  surface  and  the  electrons  relax,  due  to  the  use  of  the  localized 
basis  functions  and  the  partition  of  the  surface.  Ignoring  this  relaxation  would 
introduce  an  error.  Although  the  error  is  small  and  does  not  effect  the  dynamics 
of  electronic  states  to  a  great  extent  at  small  distance,  the  error  accumulated  at 
large  distances  could  cause  problems.  This  error  induced  by  the  partition  could 
be  significant  in  the  case  of  collisional  neutralization  of  ions. 

To  represent  the  relaxation  of  the  electrons  in  the  localized  basis  functions 
we  add  to  the  full  Hamiltonian  a  correction  term 


Vloc  =  22  (€F  ~  e™)\Xm)(Xm\ 


(6.66) 
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Figure  6.11  Pseudopotentials  for  3s  and  3p  orbitals  of  the  Na  atom. 
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where  the  summation  is  over  the  primary  region  p  and  eA  =  (xm\HM\Xm)  with 
Hm  =  —  j  V2  +  Vm  the  Hamiltonian  of  the  slab.  This  term  is  constructed  such 
that  the  asymptotic  values  of  the  electronic  Hamiltonian  are  correct.  With  this 
correction  term,  the  effective  Hamiltonian  in  the  primary  region  is 

H  =  H  +  Vloc  (6.67) 

which  will  be  used  in  the  TDHF  calculation  in  place  of  H  . 
6.5.3  Hamiltonian  Matrix  Elements 

In  the  primary  region  the  matrix  elements  of  the  Hamiltonian,  Eq.  (6.67),  are 

(6.68) 


or  explicitly, 


m£p 


H<xn  =  Han  +  (^F  ~  ^n)Sa 


(6.69) 


(6.70) 


nn  —  -"nn  T  ^F        en 


(6.71) 


"nn'    —    HsSI 


n  ^  n 


The  matrix  elements  of  H  are  calculated  as  follows 

Haa>  =  {a\-l-V2  +  VA  +  VM\a) 

=  eaSaa'  +  (a\VM\a') 


(6.72) 


(6.73) 
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Han  =  H  -  ^V2  +  VA\n)  +  (a\VM\n) 
=  eaSan  +  (a\VM\n) 
Hnn>  =  (n\  -  l-V2  +  VM\n')  +  (n\VA\n!) 


(6.74) 


(6.75) 
=  HM,Hn'  +  (n\VA\n) 

where  Em,™'  =  ("I^m|"')  does  not  depend  on  RA.  The  matrix  elements 
(a\VM\a),  (a\VM\n)  and  (n\VA\n)  depend  on  RA  and  are  calculated  by  using 
the  basis  functions  as  linear  combinations  of  Gaussians.  Neglecting  all  the  three- 
center  integrals,  we  further  have 

(n\VA\n')=Snii,{n\VA\n)  (6.76) 


When  it  approaches  the  surface,  the  ion  causes  charge  rearrangement  and 
forms  an  image  potential  in  the  solid,  which  is  important  at  low  collisional 
energies.  When  the  image  potential  is  considered  the  corrected  Hamiltonian  is 


H'  =  H  +  t>r, 


(6.77) 


where 


Vrearr  =  Vloc  +  K'r 


(6.78) 


Vim  = 


ic^-eH 


(6.79) 


CHAPTER  7 
INTEGRATION  OF  LINEARIZED  TDHF  EQUATIONS 


The  linearization  procedure  introduced  in  Ch.  4  has  been  implemented  in  a 
computation  code  and  has  been  used  in  the  integration  of  TDHF  equations  for 
both  atomic  systems  [Runge  et  al.,  90]  and  atom-surface  systems  [Feng  et  al., 
91].  In  its  applications,  we  repeatedly  tested  and  improved  this  procedure  and 
the  computing  programs.  In  this  chapter  we  discuss  some  of  the  computational 
aspects  of  our  application  of  this  procedure. 

In  Sect.  7.1  we  introduce  an  approximation  for  the  driving  term  and  present 
the  algorithm  for  integration  of  the  linearized  TDHF  equations.  In  Sect.  7.2  the 
flowchart  of  the  computation  program  for  the  integration  is  presented.  In  Sect. 
7.3  we  discuss  the  stability  and  convergence  of  the  solutions  of  the  equations  and 
the  choice  of  several  computational  parameters. 

7.1  The  Algorithm  for  Numerical  Integration  of  TDHF  Equations 


The  linearization  procedure  splits  the  TDHF  equation  into  two  linearized  time- 
dependent  differential  equations  in  a  small  time  interval  t0  <  t  <  tlt  one  for  a 
time-dependent  density  matrix  reference,  P°(/),  and  another  for  the  change  of  the 
density  matrix  due  to  the  motion  of  the  nucleus,  Q(t). 
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Omitting  the  electron-electron  interaction,  the  coefficient  matrix  W  of  the 
equation  is  time-independent  and  the  equation  for  P°(r),  Eq.  (4.37),  has  the 
formal  solution  of  the  form  of  Eq.  (4.47).  The  equation  for  Q(r),  Eq.  (4.38),  has 
a  formal  solution  of  Eq.  (4.50)  with  a  kernel  Tk^pl  which  is  a  time  integral  of 
the  driving  matrix  D'(t).  Since  the  time  interval  is  chosen  to  be  small,  we  can 
approximate  the  driving  term  by  a  function  linear  in  time,  that  is, 

D'(t)  =  C-(t-t0)  (7.1) 

where  the  coefficient  matrix  C  is  given  by 

„      D'(<i) 

C  =  —4^  (7-2) 

h  —  to 

The  term  D'(/0)  does  not  appear  in  the  above  equations,  since  the  matrix  t)'(t) 

is  proportional  to  the  changes  in  the  time  interval,  AS-1  and  AH,  and  D'(t0)  is 

actually  zero,  as  can  be  seen  from  Eq.   (5.16).   Substituting  Eq.   (7.1)  into  Eq. 

(4.51),  we  have 

FkXplit,  t0)  =  i -exp[-i(wk  -  wt)(t  -  to)] 

w/g  —  wT 

1  (7.3) 

+ -^{exp[-i(wk  -  wf)(i  -  t0)]  -  1} 

{wk  -  wj) 


7.2  Computation  Program 


All  the  computation  programs  have  been  tested  and  run  on  SUN3/50  work- 
station, SUN4/490  and  FPS  500EA1/1  computers.  The  flowchart  of  the  main 
program  of  integration  for  linearized  TDHF  equations,  TDLINEAR.F  is  given  in 
Fig.    7.1. 
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Fig.  7. 1  Flowchart  of  program  TDLINEAR.F  for  integration  of  the  linearized  TDHF  equations. 
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7.3  Stability  and  Convergence  of  the  Numerical  Integration 


The  solutions  of  TDHF  equations  are  obtained  by  integration  of  the  linearized 
equations  in  steps.  Since  possible  errors,  induced  by  the  linearization  and  the 
approximation  that  the  driving  term  is  proportional  to  time,  depend  on  the 
integration  time  interval,  integration  has  to  be  performed  in  many  small  steps. 
This  could  create  a  problem  in  that  the  computational  errors  may  be  accumulated 
in  the  course  of  the  integration  and  cause  the  solution  unstable.  The  programs 
we  used  give  stable  results  against  small  changes  of  physical  and  computational 
parameters,  such  as  atom-surface  interaction  parameters,  collision  energy,  time, 
and  initial  and  final  distances.  We  have  introduced  several  computation  parameters 
to  test  and  improve  the  stability  and  convergence  of  the  computation. 

7.3.1  Tolerances 

To  achieve  the  minimum  error  and  better  efficiency,  the  integration  is  con- 
ducted with  a  variable  time  interval  At=ti-t0,  where  r0  is  the  starting  time  and  ft 
is  the  ending  time  of  the  interval.  After  the  integration  over  one  step,  the  quotient 
HQpp(*i)||2/||Ppp(*i)||2  is  checked  against  a  designated  higher  tolerance  rh.  If 
it  is  found  greater  than  r*,  the  program  returns  to  the  original  starting  point  t0  to 
repeat  the  calculation  with  a  step  size  reduced  to  Af/2.  Such  procedure  is  repeated 
after  each  integration  until  the  quotient  is  less  than  rh,  then  the  program  proceeds 
to  perform  integration  in  the  next  step  with  the  reduced  step  size.  To  accelerate 
the  calculation  in  the  region  where  the  solutions  change  little,  we  introduce  also 
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a  lower  tolerance  rh  such  that  if  ||QPp(*i)||2/l|Pjp(*i)ll2  <  rh  the  step  size  in 
the  next  step  will  be  doubled.  Both  th  and  r\  have  been  carefully  checked  to 
assure  the  convergence  of  the  solutions. 

We  use  r/,=10^  and  r/=10-7  in  our  calculations.  Checking  calculations  with 
smaller  tolerances  7-/,=10~"7  and  r/=10~8  have  also  been  conducted  for  selected 
energies.  While  the  computing  time  for  the  smaller  tolerances  is  increased 
significantly,  the  solutions  for  the  TDHF  equations  change  little  in  the  entire 
collision.  For  example,  The  calculations  with  r/,=10~7  and  t/=10~8  required 
5,830,  3,741  and  2,156  steps  for  the  collisional  energies  of  1.0,  10.0  and  100.0  au, 
respectively;  the  calculations  with  rA=10"7  and  r^lO-8  required  11,151,  7,282 
and  4,156  steps,  respectively,  while  the  electronic  populations  changes  less  than 
0.005  au  for  all  times. 

7.3.2  Initial  and  Final  Distances 

Although  the  atom  and  the  surface  are  regarded  as  uncoupled  when  the 
distance  between  them  is  greater  than  the  interaction  range,  the  initial  and  the  final 
distances  used  in  the  calculation  must  exceed  the  distances  chosen  by  physical 
considerations. 

We  have  checked  the  convergence  of  the  results  while  changing  the  initial 
position  Zm  and  final  position  Z/  for  selected  collisional  energies,  and  found  that 
the  results  became  stable  when  Z,„  >12au  and  Z/  >16  au.  For  Z,„  and  Z/  in 
the  range  of  16-30  au,  the  electronic  populations  vary  less  than  0.005  au  for 
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collisional  energies  of  1.0,  10.0  and  100.0  au.  We  use  Zin=20  au  and  Z/=20  au 
in  our  calculations. 


CHAPTER  8 
APPLICATIONS 


The  time-dependent  molecular  method  developed  in  Chapter  2  has  been  ap- 
plied to  ion-atomic  systems  to  study  collisional  charge  transfer  and  has  provided 
useful  insight  [Runge  et  al.,  90].   In  this  Chapter  we  further  extend  our  TDMO 
method  to  ion-surface  systems  [Feng  et  al.,  91],  with  the  help  of  the  linearization 
procedure  developed  in  Chapter  4.    Ion-surface  systems  are  much  more  com- 
plicated than  atomic  systems,  and  the  direct  application  of  the  TDMO  method 
to  these  systems  will  lead  to  an  infinite  number  of  time-dependent  differential 
equations.  An  attempt  has  been  made  to  use  a  TDMO  method  in  atom-surface 
scatterings  by  modeling  the  surface  with  a  cluster  [Olson  and  Garrison,  85].  The 
partition  method  for  extended  systems  and  the  localized  basis  functions  we  have 
developed  provide  a  convenient  way  to  treat  extended  systems  and  requires  us  to 
solve  only  a  few  effective  TDHF  equations. 

The  application  presented  in  this  chapter  is  a  combination  of  the  methods 
developed  in  previous  chapters.  Our  purpose  is  to  test  the  TDMO  method  and 
partition  procedure  on  an  extended  system  and  to  develop  a  complete  and  practical 
method  for  investigation  of  dynamics  in  atom-surface  collisions.  In  order  to 
simplify  the  problem,  it  is  necessary  to  use  some  approximations  in  the  treatment 
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of  the  model,  trajectory  and  electron-electron  interaction,  which  will  be  elaborated 
in  the  text. 

In  Sect.  8.1,  we  describe  the  model  system  of  Na-W(llO)  and  calculate 
the  electronic  couplings  and  atom-surface  interaction  potentials.  Sect.  8.2  is 
devoted  to  the  trajectory  of  the  colliding  nucleus.  A  simple  potential  is  used 
in  the  application,  which  allows  an  analytic  form  for  the  trajectory.  Dynamic 
aspects  of  the  system,  such  as  initial  conditions  and  incident  parameters,  are  also 
discussed.  In  Sect.  8.3  we  describe  the  charge  transfer  process  in  the  collision. 
In  Sect.  8.4  we  present  the  results  for  the  charge  transfer.  The  evolutions  of 
the  electronic  populations  of  atomic  states  and  of  the  surface  site  are  obtained  by 
performing  the  linear  integration  of  the  TDHF  equations.  The  charge  transfer  as  a 
function  of  collision  energy  is  calculated.  To  further  explore  the  cause  underlying 
its  energy  dependence,  the  electronic  populations  are  examined  as  functions  of 
the  inverse  of  the  initial  velocity  of  the  colliding  ion. 

8.1  Na-W(llO)  Model  System 

8.1.1  Hamiltonian  and  Basis  Functions 



The  system  we  consider  consists  of  a  jellium  slab  described  in  Ch.  6  and  a 
Na  atom.  The  slab  is  located  between  -D  <  z  ^  0,  where  D  is  the  thickness 
of  the  slab.  The  location  of  the  Na  atom  is  RA(t),  which  is  a  function  of  time. 
Using  the  partition  method  developed  in  Ch.  5,  we  divide  the  system  into  a 
primary  region  and  a  secondary  region,  with  the  primary  region  consisting  of  the 
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Na  atom  and  only  one  site  on  the  surface,  which  is  the  one  at  the  "impact  point" 
at  R$  =  (0, 0, 0),  and  the  secondary  region  being  the  remainder  of  the  surface. 

The  effective  one  electron  Hamiltonian  in  the  primary  region  is  given  by  Eq. 
(6.67).  The  potential  of  the  slab  Vm  is  assumed  to  be  a  step  potential  given  by  Eq. 
(6.31)  with  a  potential  well  depth  V0=0.846  au.  A  pseudopotential  for  the  atomic 
core  is  used  for  the  atomic  potential  VA,  as  described  in  Ch.  6,  to  eliminate  the 
inner-shell  electrons.  To  simplify  the  calculation,  the  electron-electron  interaction 
is  ignored. 

The  set  of  localized  basis  functions  includes  the  generalized  Wannier  functions 
for  the  slab  constructed  in  Ch.  6  and  atomic  functions  for  3s,  3px,  3py  and  3pz 
orbitals  of  Na,  which  are  the  STO-3G  functions  described  in  Ch.  6. 
8.1.2  Electronic  Couplings 

We  neglect  the  image  potential  and  the  differential  overlap  between  the  atomic 
basis  functions  and  the  basis  functions  in  the  secondary  region,  so  that 


SaM  .  (AS"1)^  =  HaA=  AHa^=  0        m  e  s 


(8.1) 


where  s  stands  for  the  secondary  region.  The  matrix  elements  of  Hamiltonian,  Eq. 
(6.67),  are  calculated  in  Ch.  6.  The  off-diagonal  Hamiltonian  matrix  elements 
between  the  atomic  orbitals  and  between  the  atomic  orbitals  and  surface  site  R* 
are  related  to  the  electronic  couplings,  which  determine  charge  rearrangement 
during  the  collision.  In  Fig.  8.1  the  matrix  elements  HL,0  ,  H'  -  and  H'  - 
are  plotted  as  functions  of  the  distance  between  the  atomic  core  and  the  surface, 
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with  both  the  atomic  core  and  the  slab  site  of  the  primary  region  on  the  z-axis. 
The  couplings  between  3px  or  3py  and  0  are  zero  by  symmetry  and  are  not  plotted. 
At  large  distances  from  the  surface,  the  atom  behaves  just  like  an  isolated  one  and 
the  off-diagonal  elements  between  the  atomic  orbitals  of  different  symmetry  are 
zero.  As  the  atom  gets  close  to  the  surface  the  influence  of  the  surface  begins  to 
be  felt  and  the  amplitudes  of  these  elements  begin  to  increase.  The  figure  shows 
that  the  range  of  the  electronic  couplings  is  about  8  au.  The  couplings  are  large 
as  the  distance  decreases  below  4  au;  consequently,  a  strong  orbital  mixing  and 
charge  transition  may  occur. 

8.2  Atom-Surface  Interaction  Potentials  and  Trajectory 
8.2.1  Atom-Surface  Interaction  Potentials 

—  ■  — 

We  assume  that  the  atom-surface  interaction  potential  contains  contributions 
from  two  parts:  the  interaction  between  the  surface  and  the  bare  core  of  the  atom, 
and  the  average  electronic  energy.  Depending  on  the  charge  states,  the  interaction 
potential  has  different  forms.  For  the  ion  plus  surface,  the  potential  is 

Vza*  (Ra)  =  Vcc  (ilA)  +  H  (o,  o)  (8.2) 

where  Vcc  ^RA)  is  the  interaction  potential  between  the  ion  core  and  the  surface 
and  is  chosen  to  be  a  Born-Mayer  type  potential  [Hulpke  and  Mann,  83] 

Vcc(Ra)  =  Vl  Y,  exp(-a\RA  -  R*\)  (8.3) 

with  the  Born-Mayer  parameters  Vy=188.6  au  and  q=1.92  au,  where  the  summa- 
tion is  over  the  primary  region.  The  term  H'(o,  o)  =  (xff|#'|*ff)  is  the  diagonal 


104 


g 

cs 

v»-^ 

>> 

W> 

Vh 

<D 

a 

<D 

W) 

a 

•rH 

i-H 

Oh 

3 

O 

u 

1.0 


I       I       I 


-0.5 


-i 1 1 r 


i      i       i r 


-i 1 1 i '       ' 


H(3s,3pJ 

H(3s,0) 

H(3pz,0) 


■j — i — i — i — i — i ' 


0.0 


3.0  6.0 

Z(au) 


9.0 


12.0 


Figure  8.1  Coupling  energies  between  the  atomic  orbitals  and  the  surface  site  0. 
The  site  0  and  the  atomic  core  are  on  the  z-axis,  and  0  is  on  the  first  layer  at  R5. 
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matrix  element  of  the  effective  Hamiltonian  corresponding  to  the  surface  site  in 
the  primary  region.  When  an  electron  is  transferred  to  one  of  the  atomic  orbitals, 
the  atom-surface  interaction  potential  is  given  by 

Vn+A(a)  (Ra)  =  Vcc  (X)  +  H'(a,  a)  (8.4) 

where  a  is  3s,  3px,  3py  or  3pz,  and  H'(a,a)  =  {xa\H'\Xa) 

Fig  8.2  plots  the  interaction  potentials  of  V^A+  (ra\    V^^IRa), 
vY,+A(zPt)  \Ra)  and  vz+A(3Px)  (Raj  as  functions  of  the  distance  between  the  Na 
atomic  core  and  the  W(l  10)  surface.  Potential  VE+A^Py) (r\  is  not  shown  in  the 
figure  since  it  is  equal  to  VE+A^Px)  (r) .  As  one  can  see  from  the  figure,  the  con- 
tribution to  the  total  potential  from  the  electronic  energy,  which  contains  the  inter- 
action between  electrons  and  the  surface  positive  background  and  the  atomic  core, 
is  negative,  but  its  amplitude  varies  with  the  states.  In  case  of  3p*  and  3py  orbitals, 
the  electronic  energies  are  very  small  and  V^+A{3pi)(^RA^  and  VE+ A{3py)  (RaJ 
are  close  to  the  core-core  potential  Vcc(ra\     For  other  cases,  the  poten- 
tial change  due  to  the  interaction  between  electrons  and  the  slab  potential  or 
the  atomic  potential  are  strong  at  short  distance.     The  ion-surface  potential 
Vla*  \Ra)  .  which  is  below  the  atom-surface  potential  at  large  distances,  crosses 
v7:+A(33)yRAJ  and  Vfc+ Ai$p,)U&A J  &  about  z=5.5  au  and  stays  higher  as  the 
distance  decreases. 
8.2.2  Trajectory 

The  trajectory  of  the  projectile  in  the  collision  problem  has  a  significant 
effect  on  the  dynamics  of  the  electronic  states.    In  Ch.    1  we  have  discussed 
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at  length  this  subjects  as  well  as  the  merits  and  the  applications  of  a  variety  of 
approaches.  For  the  range  of  collision  energy  of  the  Na  atom  we  are  interested  in, 
the  trajectory  of  the  nucleus  can  be  treated  as  a  classical  one  which  couples  with 
the  electronic  degrees  of  freedom  of  motion.  For  example,  from  the  atom-surface 
potentials  Vj;A+yRA)  and  V^+^JrA  described  above,  one  can  construct 
average  potential  energy  curves  corresponding  to  various  configurations  to  obtain 
a  nuclear  trajectory  which  depends  on  the  charge  states  of  the  atom.  The  effect 
of  the  coupling  between  the  electronic  motion  and  the  nuclear  motion  can  also  be 
taken  into  account  by  introducing  an  effective  potential  relating  to  the  coupling 
[Runge  et  al.,  90]. 

Although  a  thorough  study  of  the  trajectory  and  its  effects  on  the  charge 
transfer  is  a  valuable  and  interesting  topic,  we  neglect  the  coupling  between  the 
nuclear  motion  and  the  electronic  degrees  of  freedom  of  motion  to  concentrate 
on  the  TDHF  method.  Instead  we  use  a  simple  trajectory  RA(t)  for  the  colliding 
atom  core  which  is  determined  by  an  atom-surface  core-core  interaction  potential 
Vec  \Ra)  ■  This  is  justified  for  the  hyperthermal  collision  energies  of  present 
interest.  For  the  same  reason  we  neglect  the  image  potential  in  the  metal. 

In  this  case  the  trajectory  has  a  prescribed  form  and  does  not  depend  on  the 
electronic  state  of  the  atom.  Let  t=0  when  ZA  is  a  minimum,  then  the  trajectory  is 


'  XA(t)  =  XA(ttn)  +  vx(t  -  tin) 
YaH)  =  YA{tin)  +  vy(t  -  tin) 


(8.5) 
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where  ma  is  the  mass  of  the  atom,  t^  is  the  initial  time,  and  vx,  vy  and  vz  are 
the  x,  y,  and  z  components  of  the  initial  velocity  of  the  nucleus,  respectively.  The 
closest  distance  from  the  surface  the  atom  can  get  to  is  decided  by  v2  and  is  given 
by  Z%tn  =  i^(^r).  In  Fig.  8.3  we  show  trajectories  of  the  Na  atom  with 
the  initial  kinetic  energies  of  1,  10  and  100  au  and  a  normal  incident  angle. 

The  initial  distance  and  the  final  distance  between  the  atom  and  the  surface 
are  chosen  such  that  further  increases  in  either  will  not  affect  charge  transfer.  It 
has  been  found  that  it  is  enough  to  use  20  au  for  both  initial  and  final  distances  for 
all  energies.  We  consider  only  the  incidence  normal  to  the  surface  in  this  study. 

8.3  Charge  Transfer 

In  the  course  of  the  collisions,  the  Na  ions  approach,  collide  with  the  surface 
and  then  recede  along  their  classical  trajectory.  In  the  meantime,  electrons  interact 
and  evolve  quantum  mechanically.  By  the  time  the  collisions  are  over,  the  charges 
have  been  considerably  rearranged,  and  some  of  them  transferred  to  the  ions.  For 
those  neutralized  ions,  a  variety  of  possible  electronic  states  may  be  populated, 
which  leads  to  the  atomic  orbital  orientation  and  alignment  For  the  normal 
incidence  by  the  ion  there  is  no  orbital  polarization  along  x-  or  y-  axes,  since 
the  only  non-zero  couplings  in  this  case  are  those  between  the  surface  state  and 
the  3s  and  3pz  orbitals. 

The  two  basic  mechanisms  for  charge  transfer  in  atom-surface  collisions  are 
the  near-resonance  process  and  the  Auger  process.  In  the  model  we  are  studying, 
since  the  energies  of  the  atomic  states  participating  in  the  charge  transfer  are  in  the 
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Figure  8.3  Trajectories  for  the  Na  atom  with  the  initial  kinetic  energies 
of  1,  10  and  100  au  and  a  normal  incident  angle.    Time  t  is  set 
to  zero  when  the  atom  is  at  the  closest  distance  from  the  surface. 
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range  of  the  electronic  band  of  the  surface  at  all  distances,  the  charge  transfer  in 
the  system  is  expected  to  be  dominated  by  the  near-resonance  process.  Thus,  we 
assume  that  the  Auger  process  is  negligible  and  only  the  near-resonance  charge 
transfer  process  takes  place  in  the  collision.  The  fact  that  all  the  potential  energies 
are  close  at  small  distances,  as  shown  in  Fig.  8.2,  justifies  our  assumption  that 
the  near-resonance  process  is  the  dominant  charge  transfer  mechanism. 

The  density  matrix  elements  for  an  isolated  surface  at  temperature  7=0K  are 


rnm  2^     UHk     km 

l*l<*F 


(8.6) 


where  the  matrix  U  is  defined  in  Ch.  6  and  kp  is  the  wave  vector  corresponding 
to  the  Fermi  energy  which  is  evalued  in  Appendix  C.  Since  the  surfaces  of  the 
jellium  slab  considered  here  are  parallel  to  the  x-y  plan,  Eqs.  (6.45)  and  (6.46) 
can  be  used  to  simplify  the  calculation. 

At  the  initial  time,  the  surface  and  the  atom  are  not  coupled  due  to  a  very 
large  initial  distance  between  them.  We  assume  both  the  surface  and  the  ion  are 
in  their  ground  states  initially;  this  leads  to  Paa'(tin)  =  0  and  PaQ(tin)  =  0  and 
the  only  non-zero  density  matrix  element  in  the  primary  region  is  P5-(£,-„)  =  P^. 

Applying  the  partition  method  developed  in  Ch.     5  and  the  linearization 
procedure  developed  in  Ch.   4,  we  come  up  with  the  effective  TDHF  equation 
(5. 15)  in  the  primary  region  with  a  linear  coefficient  matrix  Wpp.  Using  Eq.  (8.1) 
Wpp  =  (S^Ho)^  =  (S0-1)pp(Ho)pp+  (S0-1)p.(H0)lip 
=  (So-1)pP(Ho)pp 


(8.7) 


Ill 

In  the  effective  TDHF  equation  the  matrices  S,  AS,  H  and  AH  are  time  depen- 
dent through  the  known  trajectory  RaW  given  by  Eq.  (8.5).  The  equation  con- 
tains a  driving  term  D'  given  by  Eq.  (5.16),  in  which  all  the  matrix  multiplications 
are  performed  in  the  primary  region  except  for  P°8(<in)Hj8p  and  HjpsP°p(<i„), 
which  are  time-independent  and  have  to  be  calculated  only  once. 

The  Effective  TDHF  equation  is  integrated  by  use  of  the  linearization  proce- 
dure  introduced  in  Ch.  4.  An  numerical  approximation  we  used  in  this  application 
is  to  replace  the  driving  matrix  D'  by  a  function  linear  of  time  t  in  the  time  in- 
terval to  <  t  <  t\  ,  in  which  the  linearized  TDHF  equation  is  integrated.  The 
error  caused  by  this  approximation  has  been  reduced  by  decreasing  the  step  size 
until  a  stable  result  has  been  reached. 


8.4  Results 

We  have  applied  the  time-dependent  molecular  orbital  method  to  the  above 
Na-W(llO)  model  with  the  initial  kinetic  energy  of  the  atom  in  the  range  of  1 
au  to  170  au.  The  collisional  charge  transfer  is  characterized  by  the  electronic 
population  and  orbital  polarization  parameters  defined  in  Ch.  3.  For  the  normal 
incidence,  the  atomic  orbitals  are  polarized  only  along  the  z-axis. 

8.4.1  Evolution  of  Electronic  Populations 

The  evolution  of  the  electronic  density  matrix  is  obtained  by  integrating  the 
TDHF  equation  step  by  step.  The  electronic  populations  of  the  atomic  orbitals 
and  surface  site  are  calculated  according  to  the  definitions  given  in  Ch.  3.  Since 
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the  electron  spin  polarization  is  ignored,  the  electronic  population  of  an  orbital  is 
simply  twice  that  for  one  spin  on  that  orbitals. 

In  Fig.  8.4(a)  and  Fig.  8.4(b)  we  show  evolutions  of  the  electronic  populations 
of  3s  and  3pz  orbitals,  n(3s)  and  n(3pz),  respectively,  for  a  collision  energy  of 
1  au.    The  electronic  population  of  the  surface  site  0,  n(0),  is  shown  in  Fig. 
8.4(c).  In  these  figures,  the  time  is  chosen  such  that  t  is  negative  when  the  atom 
is  moving  toward  the  surface  and  is  positive  when  it  is  moving  away  from  the 
surface,  and  at  f=0  the  atom  is  at  the  closest  distance  from  the  surface.    The 
trajectory  is  plotted  in  Fig.    8.4(d)  as  a  reference  on  the  distance  between  the 
atomic  core  and  the  surface.  Initially,  the  atom  is  ionized.  At  large  distances,  the 
interaction  between  the  surface  and  the  atom  is  so  weak  that  they  behave  almost 
like  isolated  ones  without  charge  exchange  between  them.    At  around  z=5au, 
«(3s)  and  «(3pz)  begin  to  rise,  indicating  that  the  atom  enters  the  interaction 
region.    As  the  atom  gets  closer,  the  interactions  are  strong  and  the  electronic 
populations  of  the  atomic  orbitals  quickly  grow  to  significant  values.   Since  all 
the  three  orbital  energies  are  close  at  small  distances  and  even  cross  each  other 
at  some  points,  electrons  hop  among  the  three  orbitals  rapidly,  resulting  in  fast 
oscillations  in  electronic  populations.  As  the  atom  moves  away  from  the  surface, 
the  charge  exchange  is  reduced  and  the  electronic  populations  gradually  approach 
steady  values.    It  is  noticed  from  the  figures  that  the  evolution  pattern  is  not 
symmetric  on  the  way  in  and  on  the  way  out,  especially  for  3s  orbital,  reflecting 
its  dependence  of  the  charge  redistribution.  It  can  be  seen  from  the  figures  that 
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the  oscillation  frequencies  of  the  populations  change  with  the  distance.  A  careful 
check  has  shown  that  the  frequencies  are  proportional  to  the  energy  differences 
at  that  distance  but  independent  of  the  populations  of  the  states. 

Similar  patterns  are  seen  for  the  collisions  at  energies  of  10  au  and  100  au. 
Fig.  8.5(d)  shows  the  trajectory  of  the  atom  at  initial  energy  of  10  au,  Figs. 
8.5(a),  8.5(b)  and  8.5(c)  present  the  results  for  n(3s),  n(3pz)  and  n(0)  at  this 
energy,  respectively.  Figs.  8.6(a),  8.6(b),  8.6(c)  and  8.6(d)  plot  /i(3s),  /i(3pz), 
«(0)  and  the  trajectory,  respectively  for  the  collision  of  initial  energy  at  lOOau. 
8.4.2  Electronic  populations  after  collisions 

The  electronic  populations  approach  steady  final  values  as  the  atom  moves 
far  away  from  the  surface.  The  final  charge  distribution  depends  on  the  collision 
energy.  We  are  interested  in  wether  there  is  a  pattern  for  the  variation  of  the 
electronic  populations  after  the  collision  and  wether  such  a  pattern,  if  it  exists, 
can  be  understood  physically  and  provide  information  about  the  interaction  and 
dynamics  of  the  ion-surface  collision. 

We  calculate  the  electronic  populations  after  the  collision  for  a  collision  energy 
range  between  1  au  and  170  au.  The  electronic  population  of  the  3s  orbital  after 
the  collision,  «(3s),  is  shown  in  Fig.  8.7  as  a  function  of  the  collision  energy.  We 
found  that  Na  ions  are  neutralized  at  all  energy  range,  and  that  the  population  is 
relatively  lower  at  the  higher  energies,  indicating  that  the  neutralization  probability 
is  smaller.  An  energy  threshold  for  charge  transfer  of  0.02-0.04  au  has  been 
reported  for  the  Na-W  system  [Hurkmans  et  al.,  76],  below  which  the  charge 
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Figure  8.4  The  electronic  populations  and  trajectory  of  Na  atom  with  a  collision 

energy  of  1  au.  (a)  The  electronic  population  of  3s  orbital,  (b)  the  electronic 

population  of  3pz  orbital,  (c)  the  electronic  population  of  surface  site  0,  (d)  trajectory. 
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Figure  8.5  The  electronic  populations  and  trajectory  of  Na  atom  with  a  collision 

energy  of  10  au.  (a)  The  electronic  population  of  3s  orbital,  (b)  the  electronic 

population  of  3pz  orbital,  (c)  the  electronic  population  of  surface  site  0,  (d)  trajectory. 
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Figure  8.6.  The  electronic  populations  and  trajectory  of  Na  atom  with  a  collision 

energy  of  100  au.  (a)  The  electronic  population  of  3s  orbital,  (b)  the  electronic 

population  of  3pz  orbital,  (c)  the  electronic  population  of  surface  site  0,  (d)  trajectory. 
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transfer  is  cut  off.  We  found  no  threshold  in  this  study,  the  reason  being  that  our 
energy  range  is  above  the  expected  threshold.  The  electronic  population  oscillates 
between  zero  and  its  maxima  with  a  frequency  that  decreases  with  the  energy. 
There  is  an  abnormal  region  around  an  energy  of  4  au,  where  the  oscillations  are 
slower  than  that  in  the  adjacent  regions. 

In  Fig.  8.8,  the  electronic  population  of  the  3pz  orbital  after  the  collision, 
n(3pz),  is  shown  as  a  function  of  the  collision  energy.  At  small  energies,  «(3pz)  is 
zero;  it  begins  to  increase  slowly  at  E=1.9  au;  it  grows  faster  with  an  oscillatory 
amplitude  after  E=2.5  au.  A  noticed  difference  is  that  the  amplitudes  of  the 
oscillation  of  «(3pz)  are  small  although  the  frequency  is  about  the  same  as  that 
for  n(3s).  Unlike  «(3s)  whose  maxima  steadily  decreases  with  the  energy,  the 
amplitude  of  n(3pz)  varies  slowly  with  the  energy,  with  its  peak  at  4.5  au,  16  au 
and  80  au.  The  electronic  population  of  the  atom  nA,  which  is  the  sum  of  n(3s) 
and  /i(3pz)  in  our  case,  is  shown  in  Fig.   8.9. 

In  Figs.  8.10  and  8.11,  the  electronic  populations  of  3s  and  3pz  orbitals  are 
plotted  as  functions  of  the  inverse  of  the  initial  velocity,  respectively.  As  one  can 
see  from  the  figures,  the  oscillations  of  the  electronic  populations  have  almost 
constant  frequencies  in  relatively  wide  ranges  of  the  inverse  of  the  collisional 
velocity.  For  the  inverse  of  the  velocity  in  the  range  between  10  au  and  40  au, 
which  corresponds  to  energy  between  4  au  and  170  au,  the  period  of  the  electronic 
populations  is  about  1  au  of  the  inverse  velocity.  For  the  inverse  velocity  in  the 
range  between  70  au  and  150  au,  which  corresponds  to  energy  between  1  au  and  3 
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au,  the  period  of  the  electronic  populations  is  about  4.7  au  of  the  inverse  velocity. 
Finally,  in  Fig.  8.12  the  electronic  populations  of  atomic  orbitals  is  plotted  as  a 
function  of  the  inverse  of  the  initial  velocity. 
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Figure  8.7  The  electronic  population  of  3s  orbital  oi :  Na  atom 
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Figure  8.8  The  electronic  population  of  3pz  orbital  of  Na  atom 
vs.   the  collisional  energy  E.  (a)  E=l-20  au;  (b)  E=20-170  au. 
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Figure  8.9  The  electronic  population  of  Na  atom  vs.    the 
collisional  energy  E.  (a)  E=l-20  au;  (b)  E=20-170  au. 
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Figure  8.11  The  electronic  population  of  3pz  orbital  of  Na  atom  vs.  the  inverse  of 
the  collisions  velocity  v"1.    (a)  v-'=10— 35  au;  (b)  v-^35— 150  au. 
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Figure  8.12  The  electronic  population  of  Na  atom  vs.    inverse  of 
the  collisional  velocity  V1.   (a)  v'^lO— 35  au;  (b)  v-'=35— 150  au. 


CHAPTER  9 
DISCUSSION  AND  CONCLUSIONS 


This  work  presents  a  theoretical  approach  to  atom-surface  collisions  at  low 
energies.  In  Ch.  2  through  Ch.  6,  we  have  developed  a  time-dependent  molec- 
ular orbital  (TDMO)  method  to  study  charge  transfer  in  ion-surface  collisions. 
The  electronic  density  matrix  has  been  used  to  described  the  evolution  of  elec- 
tronic  states  and  distribution  of  charges  on  the  atomic  orbitals  during  and  after 
collisions.    We  have  also  used  the  density  matrix  to  define  the  orbital  polar- 
ization parameters  which  characterize  anisotropic  distribution  of  electrons  over 
orbitals  in  our  extended  systems.   To  make  the  TDMO  useful  for  applications, 
we  also  developed  other  theoretical  tools,  including  the  local  time  linearization 
method  for  integration  of  TDHF  equations,  the  partition  method  based  on  the 
localized  basis  functions  for  treatment  of  electronic  states  in  extended  systems, 
and  the  expansions  in  generalized  Wannier  functions  in  extended  systems  which 
lack  translational  symmetry.  In  Ch.  8,  we  have  applied  thie  method  to  a  model 
system.  In  this  chapter  we  discuss  these  theoretical  methods  and  their  applications 
to  the  study  of  the  dynamics  of  electronic  processes  in  atom-surface  collisions. 
Analysis  and  conclusions  are  related  to  the  physical  features  of  the  problems  and 
the  comparison  with  other  theoretical  methods.  We  also  discuss  the  application 
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of  this  approach  to  charge  transfer  in  the  collision  of  Na  with  W(110).  At  the 
end  of  this  chapter,  suggestions  for  further  work  are  furnished. 

To  study  electronic  processes  in  atom-surface  collisions,  a  pertinent  theory 
should  be  able  to  describe  both  electronic  states  and  electronic  interactions.  Some 
of  the  current  approaches  in  the  field,  while  trying  to  describe  the  evolution  of 
electronic  states,  are  unable  to  simultaneously  treat  the  electronic  interactions  in 
a  consistent  way,  and  have  to  rely  on  particular  models  or  use  parameterized 
electronic  potentials  [for  example,  Blandin  et  al.,  76].  Although  some  of  these 
approaches  come  up  with  results  comparable  with  experiments,  their  value  in 
providing  a  deep  insight  into  the  complicated  electronic  processes  seems  limited. 
The  molecular  orbital  method,  which  is  capable  of  this  task,  has  been  used 
to  study  some  time-independent  atom-surface  phenomena,  e.g.,  chemisorption, 
and  has  achieved  certain  success  [Grimley  and  Posani,  74].  For  time-dependent 
phenomena,  although  the  idea  of  the  TDMOs  has  been  used  in  some  approaches 
[Tully,  77;  McDowell,  85],  there  has  not  been  much  advancement  made  in 
obtaining  these  orbitals,  a  formal  theory  within  this  framework  has  not  been 
established  and  applied  to  collisional  phenomena  in  extended  systems. 

In  this  work,  we  have  developed  a  TDMO  method  for  atom-surface  collisions. 
The  TDMOs  are  constructed  with  atomic  functions  and  localized  basis  functions 
of  surfaces.  To  eliminate  the  artificial  couplings  at  large  distances,  introduced  by 
the  expansion  of  the  molecular  orbitals  with  atomic  orbitals  centered  on  moving 
nuclei,  the  atomic  functions  are  chosen  to  be  in  the  form  of  travelling  atomic 
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orbitals  (TAO).  The  atom-surface  and  electronic  couplings,  which  depends  on  the 
electronic  configurations  as  well  as  on  time  through  the  trajectory  of  the  atom, 
B.A(t),  are  obtained  from  the  TDMOs.  Starting  with  the  system  Hamiltonian,  we 
have  derived  the  TDHF  equation  for  the  electronic  density  matrix  which  describes 
the  evolution  of  the  TDMOs.  There  are  several  advantages  in  formulating  things 
in  this  way.    First,  by  using  the  density  matrix  we  can  correctly  describe  the 
evolution  of  the  occupied  orbitals  all  the  time  without  having  to  identify  the 
occupied  orbitals,  which  change  with  time  in  dynamic  systems.    Secondly,  the 
density  matrix  is  convenient  in  describing  the  charge  distribution  and  the  orbital 
polarization  of  subsystems.   Thirdly,  the  self-consistency  between  the  effective 
field  and  the  density  matrix  is  achieved  automatically  in  the  TDHF  equations  and 
the  self-consistent  iteration  procedure  used  in  its  time-independent  counterpart  is 
thus  avoided.  It  should  be  noted  that  the  basic  TDMO  formalism  developed  in 
Ch.  2  is  completely  general  and  will  not  rely  on  particular  basis  functions. 


One  of  the  interesting  phenomena  in  collisions  is  the  orbital  polarization 
of  the  colliding  atoms,  which  occurs  when  the  electrons  populate  the  possible 
orbitals  with  a  different  probability  and  the  electronic  cloud  has  a  distorted  form 
or  a  particular  spacial  orientation.  Study  of  the  orbital  polarization  provides 
a  deeper  and  clearer  understanding  of  the  electronic  processes  accompanying 
the  collisions.  The  TDMO  approach  can  lead  to  a  characterization  of  electron 
distribution  on  individual  states  during  the  entire  collision,  which  is  an  advantage 
over  other  approaches  which  calculate  the  total  charge  on  the  scattered  atoms  only 
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after  collisions.  We  have  defined  the  orbital  polarization  parameters  in  extended 
systems  by  tensor  operators  and  expressed  them  in  terms  of  the  density  matrix.  We 
also  defined  the  polarization  parameters  for  subsystems,  which  will  be  useful  for 
comparison  with  experimental  data  of  atomic  orbital  polarization.  As  the  distance 
between  the  atom  and  the  surface  approaches  infinity,  our  definitions  of  the  atomic 
orbital  parameters  for  the  colliding  atom  becomes  that  for  an  isolated  atom. 

Here  we  interrupt  the  discussion  of  physics  of  the  collision  problem  and  our 
method  and  turn  our  attention  to  the  mathematical  aspects  of  solving  the  TDHF 
equations.    These  differential-integral  equations,  for  which  the  exact  solutions 
are  hard  to  obtain,  are  often  solved  by  linearization  procedures,  which  linearize 
the  equations  by  assuming  that  the  changes  of  their  solutions  are  slow  in  a 
small  time  interval.  However,  since  the  density  matrix  is  expected  to  be  a  fast 
oscillatory  function  in  our  case,  we  can  not  use  the  conventional  procedure  that 
uses  a  constant  zeroth  order  solution.    We  notice  that  the  time  dependence  of 
the  density  matrix  results  from  two  causes:    the  change  of  the  effective  field 
due  to  the  nuclear  motion,  and  the  change  due  to  electronic  transition  between 
states.  While  the  former  changes  slowly  for  low  energy  collisions,  the  electrons 
jump  between  states  with  a  frequency  proportional  to  the  inverse  of  the  energy 
difference  of  the  states,  that  could  be  high  in  case  of  the  near-resonance  charge 
transfer  process.   The  linearization  procedure  we  proposed  treats  the  change  of 
the  density  matrix  in  a  short  time  as  a  sum  of  a  time-dependent  reference  density 
matrix  P°(t),  which  satisfies  a  linearized  TDHF  equation  for  the  density  matrix 
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when  the  nuclear  position  is  fixed  in  space,  and  a  matrix  Q(t),  which  is  the  first 
order  change  of  the  density  matrix  and  satisfies  another  linearized  TDHF  equation 
which  is  coupled  with  that  for  P°(<).  We  have  obtained  for  a  general  case  the 
formal  solutions  for  P°(t)  and  Q(t)  from  an  exponential  transformation,  and 
the  analytical  solutions  for  a  special  case  where  the  electron-electron  interaction 
is  ignored.    Computational  programs  have  been  developed  to  solve  the  TDHF 
equations  in  atom-surface  collisions  using  this  linearization  procedure.  The  key 
to  implementing  this  procedure  is  that  Q(t)  should  be  small  compared  with  P°(t). 
We  employ  variable  step  sizes  and  check  this  criterion  after  every  integration.  The 
stability  and  convergence  of  the  solutions  have  also  been  checked. 


One  of  the  major  problems  when  applying  the  TDMO  approach  to  extended 
systems  is  the  treatment  of  the  great  number  of  electronic  states.  Collisional 
charge  transfer  in  atom-surfaces  interactions  is  a  process  which  involves  a  few 
electrons  and  strong  couplings  in  a  short  period  of  time.  We  have  developed  a 
partition  method  to  split  the  system  into  a  primary  region  and  a  secondary  region 
and  constructed  in  the  small  primary  region  the  effective  TDHF  equations  which 
are  coupled  with  the  secondary  region  through  a  time-dependent  driving  term.  In 
the  secondary  region,  the  electrons  are  effected  by  the  colliding  atom  indirectly 
through  their  coupling  with  the  electrons  in  the  primary  region,  and  the  states  of 
the  electrons  are  approximately  unchanged  during  the  collision.  Apparently,  the 
primary  region  should  be  chosen  to  have  a  size  of  the  order  of  the  interaction 
range;  we  will  come  back  to  it  in  the  discussion  of  the  results.   This  partition 
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method  can  also  be  applied  to  other  systems  which  contains  strong  and  local 
interactions. 

Our  partition  method  requires  a  set  of  localized  basis  functions.   We  have 
further  extended  the  idea  of  using  the  generalized  Wannier  functions  (GWFs) 
for  extended  systems,  which  are  orthonormal  and  localized  around  the  sites 
in  the  systems.    We  have  derived  a  formalism  to  obtain  these  functions  by  a 
variational  principle,  which  shows  that  GWFs  can  be  constructed  from  the  system 
Hamiltonian  without  having  the  knowledge  of  the  eigenfunctions.  These  GWFs 
are  written  in  the  form  of  linear  combinations  of  Gaussians  for  practical  reasons 
and  can  be  used  in  studies  of  local  phenomena  in  extended  systems  which  lack 
translation  symmetry.  We  also  argue  that,  because  the  GWFs  quickly  recover  the 
form  of  the  Wannier  functions  of  a  bulk  when  they  are  away  from  the  defects, 
only  a  small  number  of  functions  need  to  be  determined  for  quasi-periodic  systems 
which  lose  translational  symmetry  only  in  certain  directions  or  contains  defects  in 
local  areas.  The  GWFs  have  been  constructed  for  a  jellium  surface,  for  which  the 
three  dimensional  calculation  is  reduced  to  a  one  dimensional  calculation.   The 
eigenvalues  and  eigenfunctions  calculated  using  GWFs  as  basis  functions  are  in 
good  agreement  with  the  exact  ones.  To  complete  our  basis  set,  we  constructed 
the  atomic  basis  functions  for  the  valence  orbitals  through  a  pseudopotential  for 
the  atomic  core.  In  doing  so,  the  inner  electrons  are  eliminated  since  they  do  not 
play  an  important  role  in  charge  transfer. 

We  have  applied  the  TDMO  method  to  the  collisional  charge  transfer  in 
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Na-W(l  10).  The  surface  is  modeled  by  a  jellium  with  its  band  parameters  fitted  to 
those  of  the  W(110)  conduction  band.  In  the  present  applications  of  this  method, 
some  simplifications  are  used,  including  ignoring  the  electron-electron  interaction, 
using  a  prescribed  trajectory,  and  using  a  simple  model  for  the  surface.  Discussion 
of  the  results  includes  the  electronic  couplings  and  atom-surface  interactions, 
evolution  of  the  electron  populations  and  the  energy  dependence  of  the  electron 
populations. 

We  have  calculated  the  electronic  couplings  and  the  interaction  potentials  of 
the  system.  The  couplings  between  the  atomic  orbitals  show  an  interaction  range 
of  8  au.   At  the  short  distances,  the  couplings  between  the  surface  and  atomic 
orbitals  are  about  the  same  order  as  that  calculated  from  the  WKB  approach 
in  parabolic  coordinates  [Grozdanov  and  Janve,  78]  and  that  from  muffin-tin 
potentials  [Muscat  and  Newns,  79];  there  are  no  other  results  to  compare  with 
for  the  coupling  between  the  3s  and  3pz  orbitals.     The  strong  couplings  at 
short  distances  indicate  that  the  orbitals  are  heavily  mixed  and  a  strong  electron 
exchange  is  possible.   At  larger  distances,  the  decay  of  the  couplings  is  faster 
than  the  results  calculated  by  other  methods,  which  could  be  attributed  to  the 
fact  that  only  a  small  number  of  Gaussians  are  used  in  our  calculations  when 
writing  the  basis  functions  as  linear  combinations  of  Gaussians.  The  results  for 
interaction  potentials  corresponding  to  various  configurations  are  shown  in  Fig. 
8.2.    However,  there  is  to  our  knowledge  no  experimental  data  for  Na-W  to 
compare  with.  All  the  potentials  except  for  Na(3px)-W+(0)  and  Na(3py>-W+(0) 
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have  minima  at  certain  distances  from  the  surface.  The  positions  of  the  minima 

for  Na(3s)-W+(0)  and  Na(3pz)-W+(0)  are  at  4.0  au,  while  that  for  Na+-W(0)  is 

4.6  au.  The  important  information  obtained  from  the  calculation  is  that  potentials 

corresponding  to  different  configurations  are  within  0.1  au  at  short  distances, 

which  makes  the  near-resonance  charge  transfer  process  the  dominant  one. 


The  evolution  of  the  electronic  populations  for  atomic  orbitals  and  surface 
orbital  have  been  calculated,  which  offers  a  useful  insight  into  the  dynamic 
behavior  of  charge  transfer.  Unfortunately,  spectroscopy  techniques  have  not  yet 
been  developed  to  provide  corresponding  data  for  comparison.  Common  to  all  the 
energies,  the  charges  begin  to  slowly  transfer  to  the  atomic  orbitals  when  the  atom 
approaches  to  8  au  from  the  surface,  and  then  rapidly  increase  at  about  5  au  where 
the  potential  curves  cross.  All  the  electron  populations  contain  fast  oscillations 
within  the  distance  range  of  5  au,  which  is  characteristic  of  the  near-resonance 
charge  transfer.  Recalling  that  the  density  matrix  is  the  sum  of  a  time-dependent 
background  P°(t)  and  a  small  change  Q(t),  the  rapid  oscillations  are  decided 
by  the  evolution  of  P°(t)  and  slow  changes  of  the  amplitudes  are  attributed  to 
the  accumulation  of  the  change  of  the  density  matrix  Q(t).  The  pattern  of  the 
evolution  of  the  populations  varies  with  the  orbital  and  energy.  Noticing  that  these 
curves  on  the  outward  path  have  different  shapes  than  those  on  the  inward  path, 
it  means  that  the  evolution  of  the  electron  populations  depends  on  the  amount  of 
charge  on  the  orbitals.  The  frequency  of  the  oscillations  varies  with  the  distance 
and  is  higher  in  the  short  distance  range.  Our  check  indicates  that,  for  a  collisional 
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energy,  the  frequency  is  roughly  inversely  proportional  to  the  energy  differences 
between  the  potential  energies  of  atomic  state  and  the  surface  state.  This  seems 
to  support  the  explanation  that  the  transition  rate  between  two  orbitals  depends 
on  the  energy  difference.  However,  the  frequency  also  varies  with  the  collisional 
energy,  and  is  higher  at  lower  energies,  probably  because  the  coupling  between 
the  electronic  motion  and  the  nuclear  motion  should  be  taken  into  account  at  lower 
energies.  The  total  charge  in  the  primary  region  before  and  after  the  collision  but 
has  been  found  to  change  less  than  15%  for  E=l-120  au,  which  is  an  indication 
that  the  charge  transfer  occurs  indeed  mainly  in  the  primary  region.  A  more 
accurate  calculation  should  however  consider  a  larger  primary  region,  especially 
for  high  energies. 

We  have  studied  the  energy  dependence  of  the  electron  populations  and 
found  oscillations  in  the  whole  energy  range  between  1  au  and  170  au,  with 
a  frequency  that  decreases  steadily  with  increasing  the  collision  energy.  Such 
an  oscillatory  structure  has  been  encountered  in  some  experiments  on  ion-metal 
surface  collisions  [Erickson  and  Smith,  75;  Taglauer  and  Heiland,  76;  Tolk  et  al., 
76].  This  is  believed  to  arise  from  the  quantum  phase  interference  between  the 
states  whose  energies  are  close  in  the  near-resonance  charge  transfer  process,  and 
will  be  explained  here  by  a  simple  model  [Tolk,  et  al.,  76].  Consider  an  atom- 
surface  system  with  two  potential  curves  corresponding  to  the  ion-surface  and 
atom-surface  configurations,  which  are  split  by  a  small  difference  AE(R),  where 
R  is  the  distance  of  the  atom  to  the  surface.  When  the  atom  enters  the  interaction 
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region,  the  coupling  between  the  two  states  becomes  comparable  with  AE(R)  and 
the  charge  transfer  occurs.  As  the  electrons  on  these  two  states  evolve  along  their 
potential  curves,  a  differential  phase  A$  is  developed  between  the  two  states, 
until  the  atom  leaves  the  interaction  range.  After  the  collision  the  populations  of 
the  atomic  configuration  and  ion  configuration  are  expressed  as 


na  =  aa  +  0cos2  (A<£/2) 


(9.1) 


m  =  a,  +  pcos2(A<t>/2)  (9.2) 

where  aa  ai  and  f3  are  slowly  varying  functions  of  the  incident  ion  energy.  The 
phase  is  given  by  (using  h=l) 


*->*>*-/35p* 


(9.3) 


where  v(R)  is  the  normal  component  of  the  nuclear  velocity.  Since  the  velocity 
is  nearly  constant  in  a  wide  range,  the  above  equation  is  approximated  by 


■Km 

M  «  -  /  AE(R)dR 

Vi  J 

Ra 


(9.4) 


where  u,  is  the  initial  velocity  component  perpendicular  to  the  surface,  and  Rm  and 
Ro  are  the  interaction  range  and  turning  point,  respectively.  This  Model  predicts  an 


oscillation  frequency  proportional  to  the  inverse  of  the  incident  velocity,  v'1.  We 
have  studied  the  dependence  of  populations  on  v'1  and  found  that  the  oscillation 
frequency  is  roughly  constant  in  the  region  u"1  =70- 144  au  and  in  the  region 
v~  =10-40  au.  The  fact  that  the  frequency  changes  slowly  with  v{  indicates  that 
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v(R)  is  not  constant.  To  further  check  this  theoretical  explanation,  the  calculation 
should  be  done  for  the  collisions  with  different  incident  angles  to  examine  the 
angular  dependence  of  the  peak  positions.  In  ionization  spectroscopy  experiments 
[Overbosch  et  al.,  80]  the  oscillation  of  the  ion  yield  is  not  seen.  A  possible 
reason  is  that  the  experiment  did  not  measure  enough  points.  It  is  also  possible 
that,  since  the  oscillations  of  the  populations  are  very  fast  in  the  low  energy  range, 
the  final  yield  of  the  ions  in  the  experiment  is  actually  an  average  of  the  product 
yields  with  different  initial  velocities. 

The  charge  transfer  exhibits  a  difference  between  the  3s  and  3pz  orbitals.  For 
3s,  the  charge  transfer  occurs  in  the  whole  energy  range  between  1  au  and  170  au 
and  its  intensity  oscillates  between  0  and  peaks.  In  contrast,  the  oscillations  for 
3pz  are  smaller  in  the  low  energy  range,  and  the  amplitude  of  the  peaks  is  slowly 
varying  with  three  maxima  at  E=4.5  au,  16  au  and  80  au.  The  threshold  for  Na-W 
reported  in  experiment  [Hurkmans  et  al.,  76]  is  0.02-0.04  au.   Our  calculation 
shows  that  the  3pz  orbital  has  a  very  high  threshold  at  2.5  au.    However,  the 
evolution  calculation  for  collision  energy  of  1  au  shows  that  charges  are  actually 
transferred  to  3pz  orbitals  when  the  atom  is  in  the  interaction  range.  They  jump 
between  3pz  and  other  orbitals  and  then,  on  the  outward  path,  transfer  to  the 
surface  and  3s  orbitals.  The  calculation  should  be  extended  to  lower  the  energy 
range  to  see  wether  a  threshold  for  3s  orbitals  exists.  It  will  also  be  interesting  to 
do  calculation  for  different  incident  angles  to  study  the  time-dependence  and  the 
thresholds  of  the  electronic  populations  for  3px  and  3py  orbitals. 


136 
This  work  has  constructed  the  basic  framework  of  the  TDMO  method  and 
introduced  some  accessory  tools  for  applying  it  to  atom-surface  collisions.  In 
the  development  of  this  method,  we  have  not  tempted  to  use  empirical  data  to 
fit  the  our  results  with  experimental  ones.  We  believe  that  keeping  the  theory  in 
its  most  general  form  is  probably  the  biggest  intrinsic  advantage  of  the  TDMO 
method  over  the  approaches  which  can  not  treat  the  effective  potentials  and  the 
charge  transfer  self-consistently  and  have  to  rely  on  the  empirical  results.  In  our 
applications,  we  have  made  some  approximations  and  simplifications  which  can 
either  be  removed  or  be  made  more  accurate. 

However,  insofar  our  model  is  simple  and  the  results  are  preliminary,  some 
extension  and  applications  of  this  method  are  suggested  here. 

(a)  Size  of  the  primary  region.  In  principle,  the  size  of  the  primary  region 
should  be  of  the  order  of  the  interaction  range,  but  the  experiences  of  other  works 
using  similar  partition  methods  have  shown  that  a  smaller  primary  region  could 
give  satisfactory  results  [McDowell,  82;  de  Malo  et  al.,  87].  Investigation  of  the 
convergency  of  results  with  the  primary  region  size  can  determine  the  optimum 
size  of  the  primary  region. 

(b)  Orbital  polarization.  We  have  formulated  the  method  to  study  the  orbital 
polarization.  However,  the  orbital  polarization  shows  in  a  limited  way  in  the  appli- 
cation presented  in  this  work  when  the  collision  is  limited  to  the  normal  incidence. 
Application  to  incidence  with  a  finite  angle  will  provide  more  information  about 
the  electronic  interactions  in  the  collisions,  and  will  be  a  new  test  of  the  theory. 
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(c)  Trajectory.  Again,  due  to  the  limit  of  time,  the  application  of  the  method 

has  been  confined  to  use  a  simple  prescribed  trajectory.  Although  the  trajectory 
study  is  relatively  independent  of  the  TDMO  method  itself,  it  should  be  given  a 
great  deal  of  attention  for  its  importance  to  the  problem.  Our  studies  should  be 
extended  to  include  the  coupling  between  the  electronic  motion  and  the  nuclear 
motion.  Such  a  study  could  use  the  effective  force  containing  the  couplings,  along 
the  lines  we  have  applied  to  atomic  collisions  [Runge  et  al,  91],  or  within  the 
eikonal  method  [Micha,  83]. 

(d)  Electron-electron  interaction.  Ignoring  the  electron-electron  interaction 
in  our  application  is  not  necessary,  but  it  simplifies  the  calculations.  Including 
this  interaction  will  make  the  TDHF  equations  more  demanding  of  computing 
time.  Another  aspect  of  electron-electron  interactions  is  the  appearance  of  image 
potentials  in  the  metal.  These  can  be  neglected  at  the  higher  energies  we  consider, 
but  would  be  more  important  in  threshold  studies. 

(e)  Basis  functions.  The  TDMO  method  has  been  cast  in  such  a  way  that  it 
uses  localized  basis  functions  but  does  not  depend  on  the  particular  form  of  them. 
This  allows  the  freedom  of  applying  this  method  to  more  practical  models  and 
using  other  localized  basis  functions.  An  example  is  the  localized  basis  functions 
developed  by  Smith  and  his  co-workers  in  studies  of  surface  electronic  structure 
[Smith  and  Gay,  75;  Smith  et  al.,  80]. 


APPENDIX  A 
CALCULATION  OF  COEFFICIENT  MATRIX  Bx 


The  one  dimensional  Wannier  function  for  a  free  electron  gas  is  [Callaway, 


76] 


1       sm[%(x  -  nsd)] 
wnx{x)  =  ~7j FT" n     » 


(A.1) 


where  d  is  the  lattice  constant.  We  consider  a  function  Wn'xm'(x)  which  is  a  linear 
combination  of  2LX+1  Gaussians 


mx=nx+Lx 


(A.2) 

where 

jj,(^„i)  =  &m,exp  -/?TO,(x  -nrd) 

is  a  one  dimensional  Is  primitive  Gaussian  with  a  constant  bmz  =  ( ^u-) 

We  use  the  least  square  method  [Shavitt,  63]  to  determine  the  coefficient  ma- 
trix B%xmx(0mz)  which  minimizes  the  difference  between  w„x(x)  and  Wn™*^). 
The  deviation  between  w*  (x)  and  Wn™'{x)  is 


(A.3) 

1/4 


to 


mx=nx+Lx 


dx 


(A.4) 


The  choice  of  the  linear  coefficients  B%am  (/9m.)  which  leads  to  the  minimum 
D,  according  to  the  least  square  method,  requires  the  vanishing  of  the  partial 
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derivatives,  dD/dB£xTnx  =  0,  which  gives 

9mx(Pmz,x) 


I 


m'x=nx+Lx 


<fx  =  0         (A.5) 


Defining  the  overlap  matrix  of  the  Gaussians 


GSU.(A»J  -  WL(Afc)WL0^)) 


(A.6) 


and  a  matrix 


<£,m,(fl».)  -  «fo£.(/0) 


(A.7) 


Eq.   (A.5)  is  then  written  as  a  matrix  equation 


Cx  =  BxG, 


(A.8) 


from  which  we  obtain  the  coefficient  matrix 


X//~lX\  —  1 


Bx  =  CX(GX) 


(A.9) 


Because  of  the  translation  symmetry  in  the  x  direction 


Bnxmx(Pmt)  =  £|;OT#(A»J 


and  only  one  row  of  Bx  must  be  found. 


(A.10) 
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where  fim  =  £$£-,  Rm  = 


"S        "m 


and 


(B.2) 


APPENDIX  B 
CALCULATION  OF  TOTAL  BAND  ENERGY  ft 


For  a  jellium  slab  with  a  step  potential 

'0  z>0 ' 

K(f)  =  I  -Vo  -D<z<0    \  (B.l) 

.0  z<-Dd 

the  total  energy  of  a  band,  using  the  generalized  Wannier  functions  u;^(f)  which 

are  written  as  of  linear  combinations  of  Gaussians  g^(/3mz,f),  is 

* 

=  E  E  E  B*aB*f  >M$mM\M>{M) 

n       m     m' 

where 

is  an  element  of  the  coefficient  matrix  B,  and  H  is  the  Hamiltonian, 

1  n 

H  =  -3  V2  +  V(f) 
The  matrix  element  for  the  kinetic  energy  is  [Shavitt,  63] 

(<7nl  "  g  V2|^m)  =  fe(3  -  Wn*Rlrn)Gnm 


(B.3) 


(B.4) 


(B.5) 


Gm  =  {9n\gA)  =  (^4^m  )     ^p(-fe^L)  (B.6) 
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is  the  overlap  matrix  of  the  Gaussian  functions. 

Noticing  the  potential  V(r)  is  constant  in  x  and  y  directions,  the  matrix 
elements  for  the  potential  are 

(9n\V\9rn)  =  (rf,(A».M,0?m,))  ■  <<(/?mJ| < (/?m J) 

•b..\V\gm.)  (B.7) 

=  QLviflm.)  ■  OUttfU,)  •  {9nM\V\gm.) 


Calculation  shows  that 


Vt 


(9n,\V\gmz)  =  -^-[erf(l2)-erf(l1)} 

where  erflj)  is  the  error  function  [Gauatschi,  72]  and 

'/3nxnz  +  f3mzmz 


h  =  VPn,  +  0m, 


0n,  +  Pm, 


(B.8) 


(B.9) 


H  =  VPn,  +  Sb. 


,0  + 


j3n,nz  +  ftro,ma 
At.  +  An. 


d 


(B.10) 


Substituting  Eqs.  (B.5)  and  (B.  8)  into  Eq.  (B.4)  the  matrix  element  of  the 
Hamiltonian  becomes 

(9*\  ~  |v2ifta)  =  $m{Am()  -  qfaAJD 


-  y[er/(/2)  -  «•/(/!)]} 


(B.ll) 


Using  Bx  and  By  obtained  in  Appendix  A  and  substituting  Bz  =  (Gz)    2  and 
Eq.  (B.ll)  into  Eq.  (B.l),  the  total  band  energy  Q  can  be  calculated. 
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APPENDIX  C 
CALCULATION  OF  ELECTRONIC  PROPERTIES  OF  THE  FINITE  SLAB 


We  consider  a  finite  jellium  slab  which  models  the  W(110)  surface  and  has 
Ns  =  NxNyNz  sites  which  form  a  cubic  mesh,  where  Nx  is  the  number  of  sites  in 
the  x-direction,  Ny  is  the  number  of  sites  in  the  y-direction  and  Nz  is  the  number 
of  sites  in  the  z-direction.  We  choose  Nx=21,  Ny=2l  and  Nz=22.  It  should  be 
pointed  out  that  the  sites  are  chosen  from  mathematical  considerations  and  are  not 
the  lattice  points  in  W(110)  and  that  the  distance  between  sites,  d,  is  not  equal 
to  Tungsten's  lattice  constant  a.  In  Sect.  6.2  we  calculated  generalized  Wannier 
functions  for  the  slab,  assuming  it  satisfies  the  cyclic  boundary  conditions  in  x- 
and  y-directions.  In  this  appendix  we  describe  how  the  Fermi  energy  of  the  slab 
was  calculated. 

Tungsten  has  an  FCC  structure  with  a  lattice  constant  a=5.97  au;  in  each  unit 
cell  there  are  two  atoms  with  six  valence  electrons  each,  which  gives  a  density 
of  valence  electrons  p=0.0564  au.  We  choose  d=2.6l  au  so  that  average  number 
of  electrons  per  mesh  cell,  or  pd=l. 


The  energy  of  a  level  is 


E%  =  Ekx  +  Eky  +  Ekz 


(CI) 
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with  Ekx  =  &,  Eky  =  f ,  where 


tCx  — 


2njx 


(Nx  -  l)d 


jx=  0,1, 


(C.2) 


_       2xjy 
ky  ~  (N,  -  l)d 


jy=0,l, 


(C.3) 


and  Ek,  is  the  exact  eigenvalue  for  kz(jz),  jz=l,  2,  .  .  .,  for  a  one-dimensional 
square  well,  which  is  listed  in  Table  6.2. 

Starting  with  the  lowest  level,  we  fill  each  level  with  two  electrons.  When 
all  the  N=(Nx-l)(Ny-l)(Nz-l)  electrons  are  located  in  levels,  the  energy  of  the 
highest  occupied  level  is  the  calculated  Fermi  energy.  Using  the  bottom  of  the 
band  eb=-0-846  au,  we  found  that  the  calculated  Fermi  energy  is  eF=^-0-192  au, 
which  is  in  agreement  with  the  experimental  value  of  W(110),  -0.191  au. 
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